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ABSTRACT 

We have developed a new computer code, RAM, to solve the conservative equations of special rel- 
ativistic hydrodynamics (SRHD) using adaptive mesh refinement (AMR) on parallel computers. We 
have implemented a characteristic- wise, finite difference, weighted essentially non-oscillatory (WENO) 
scheme using the full characteristic decomposition of the SRHD equations to achieve fifth-order ac- 
curacy in space. For time integration we use the method of lines with a third-order total variation 
diminishing (TVD) Runge-Kutta scheme. We have also implemented fourth and fifth order Rungc- 
Kutta time integration schemes for comparison. The implementation of AMR and parallelization is 
based on the FLASH code. RAM is modular and includes the capability to easily swap hydrodynamics 
solvers, reconstruction methods and physics modules. In addition to WENO we have implemented a 
finite volume module with the piecewise parabolic method (PPM) for reconstruction and the modified 
Marquina approximate Ricmann solver to work with TVD Runge-Kutta time integration. We exam- 
ine the difficulty of accurately simulating shear flows in numerical relativistic hydrodynamics codes. 
We show that under-resolved simulations of simple test problems with transverse velocity compo- 
nents produce incorrect results and demonstrate the ability of RAM to correctly solve these problems. 
RAM has been tested in one, two and three dimensions and in Cartesian, cylindrical and spherical 
coordinates. We have demonstrated fifth-order accuracy for WENO in one and two dimensions and 
performed detailed comparison with other schemes for which we show significantly lower convergence 
rates. Extensive testing is presented demonstrating the ability of RAM to address challenging open 
questions in relativistic astrophysics. 

Subject headings: hydrodynamics - methods: numerical - relativity 



1. INTRODUCTION 

Many astrophysical phenomena involve gas moving at 
relativistic speeds. Classical sources include active galac- 
tic nuclei (AGN), microquasars, pulsar wind nebulae and 
gamma-ray bursts (GRBs). More recently, conclusive ev- 
idence has mounted that a subset of core collapse super- 
novae !aje_assj3ci^ed^tli^^ 5 s) GRBs 
re.g..lHiorth et~aTll2003lSta,nek et al]l20T)l. Astrophys- 
ical processes at the endpoint of stellar evolution thus 
appear capable of accelerating flows to ultra-relativistic 
speed (W > 100, where W denotes Lorentz factor.) 

Additionally, the fading afterglows of cosmological 
GRBs are now o bserved with an incr easing rate by the 
SWIFT satellite <|Gehrels et al]l2004l) . GRB afterglows 
are produced after the gamma-ray producing relativis- 
tic flow transfers its energy to the circum-burst medium 
in the form of a strong relativistic shock. Of particular 
interest is evidence that the ejecta producing GRBs and 
their afterglows are beamed into jets. The quality of cur- 
rent afterglow observations require high-resolution multi- 
dimensional simulations of jetted relativistic shocks for 
interpretation. 

Relativistic hydrodynamics simul ations are more dif- 
ficult than Newtonian simulations ijNorman fc W inkler 
1986). However, there has been much progress in rela- 
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tivistic numerical simulations during the past decade due 
to the development of modern numerical methods and 
the increasing power of computers. We refer the reader 
to a comprehens ive review of numerica l relativistic hy- 
drodynamics by iMarti fc Miillerl l)1999t) and references 
therein. The so-called high-resolution shock-capturing 
(HRSC) methods, which are based on the fact that 
special relativistic hydrodynamics (SRHD) with causal 
equat ion of state is a hyperbolic system of conservation 
laws (j Anild 1 1 989|) . are particularly promising in modern 
numerical simulations. They can achieve very high ac- 
curacy and handle ultra-relativistic flows, strong shocks 
and contact disco ntinuities extremely well. GENESIS 
(lAlov et alJ IT^ is a very efficient scheme based on 
HRSC techniques. 

The high-order essen tially non-oscillato ry (ENO) 
shock capturing schemes l|Harten et al.lll987l) have been 
very successful in numerically solving hyperbolic systems 
e.g., the Eulerian gas dynamics equations. ENO-based 
methods have previously been im plemented for computa- 
tional relativistic hydrodynamics dDolezal fc Won ell995t 
iDonat et allll998l: iDel Zanna fc Bucciantinill2002|) . The 
weigh ted essentially non- oscilla tory (WENO ) sch emes 
Ce.g.. ILiu. Osher fc ChanllT99l Ijiang fc Shul 119961) are 
recent developments following the philosophy of ENO 
schemes. 

Adaptive Mesh Refinement (AMR) has played an in- 
creasingly important role in many branches of numerical 
astrophysics including e.g., cosmology, star formation, 
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jets, interacting binaries, stellar wind collisions , and su- 
perno va explosion and remnant evolution (see iNormanl 
12004 and references therein). Examples of recent work 
in numerical astrophysics inclu de using AMR to solv e 
the equations of hydrodynamics (iPlgwa&M^illcr 2001), 
magnetohydrodynamics (MHD) (Balsara il200lD. specia l 
relativistic hydrodynamics (SRHD) (Hughes et al. 2002), 
general relativistic hydro dynamics (|Donmez| I 2004D and 
general relativistic MHD (Anni nos et al]l2005j) . 

In this paper, we describe a new modular, highly 
accurate, special relativistic hydrodynamics code with 
adaptive mesh refinement, RAM (Relativistic Adaptive 
Mesh). The modular design of the code allows us to 
easily change between various algorithms. We have im- 
pleme nted the fifth-order WENO scheme of lJiang k Shul 
( 1996) in SRHD for the first time as one of the modules. 
We have also implemented an H RSC module using the 
modified Marqu ina flux formula ijMarauina et al.lll992t 
| Alov_e^an[1999j) and piecew ise parabolic reconstruction 
llColella k Woodward! IT984) similar to the GENESIS 
code ijAlov et al.lll999|) . We also use piecewise linear re- 
construction (PLM) with both WENO and HRSC and 
perform comparisons between all methods for several 
standard test cases. We present detailed tests of the code 
demonstrating its ability to handle the extreme resolu- 
tion requirements of simulating ultra- relativistic flows. 

A primary motivation in writing the RAM code is to 
study the relativistic explosions producing cosmological 
GRBs. We plan to use RAM to simulate the relativistic 
flows produced at the endpoint of stellar evolution. In 
particular, numerical simulations of the collapsar mode l 
for GRBs l)Wooslevlll99l iMacFadven k WooslevllT999]) 
are very challenging because of the ultra-relativistic 
speed and detailed micro-physics involved in the prob- 
lem. Relativistic hydrodynamical simulations of jets in 
collap sars have been done successfully by sever al au- 
thors ijAlov et al] 120001 iZhang et all 1200.4 12004) using 
the GENESIS method. We plan to use RAM to extend 
these calculations to higher resolution over a larger dy- 
namic range in lengthscale made possible by AMR. 

In §Elwe present the fundamental equations of SRHD 
and describe the conserved variables evolved by our code. 
In § |3 we describe the flux differenced semi-discrete 
form of the equations we solve, the algorithms used for 
construction of numerical fluxes and how the equations 
are integrated in time. In § 0] we present results from 
standard test problems on a fixed uniform mesh includ- 
ing convergence tests. In § we describe the adap- 
tive mesh refinement algorithm employed in RAM by 
utiliz ing components of the FLASH code ijFrvxell et alJ 
120001). which in turn a dapted the PARAMESH package 
ijMacNeice et al.ll2000l) . In § [S] we present test problems 
run with AMR in one, two, and three dimensions includ- 
ing Riemann problems with transverse velocity, a 2D ax- 
isymmetric jet in cylindrical coordinates and a 3D blast 
wave. In §0we summarize the results of the paper. In 
Appendix ^ we provide details of our implementation of 
SRHD in curvilinear coordinates. 



2. GOVERNING EQUATIONS OF SPECIAL 
RELATIVISTIC HYDRODYNAMICS 

The governing equations of special relativistic hydro- 
dynamics (SRHD) describe the conservation of rest mass 



and stress-energy of a fluid: 

(p^);» = 0, (1) 

and 

(T^)-m = (2) 

where p is the rest mass density measured in the fluid 
frame, u M = W(c, u) is the fluid four-velocity, W is the 
Lorentz factor, c is the speed of light, u is the classical 
three- velocity, T^ v is the stress-energy tensor of the fluid 
and the subscript ;M denotes the covariant derivative. For 
a perfect fluid the stress-energy tensor is 



T» v = phu^u" + pg 



(3) 



where h = 1 + e + p/pis the relativistic specific enthalpy, 
e is the specific internal energy, p is the pressure and 
g^ v is the inverse metric which here we take to be the 
Minkowski metric though extension to curved spacetimes 
is evident. 

SRH D can be written as a set of conservation laws (see, 
e.g. IMarti k MuIIctIITqQQT) , 
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where the conserved variable U is given by 

U = (D,S\S 2 ,S 3 7 t) t , (5) 
and the fluxes are given by 

F J = (Dv j ,S 1 vi+p6 j 1 ,S 2 v j +p6 j 2 ,S\ j +p6 j 3 ,S : >-Dv j ) T , 

(6) 

here (5 M „ is the Kronecker symbol and v° is the velocity. 
The conserved variables U include rest mass density, D, 
the momentum density, S 3 , and energy density, r. They 
are measured in the laboratory frame, and are given by 
(assuming the speed of light c = 1), 



D = pW 
S j =phW 2 v 3 
r = phW 2 -p- pW, 



(7) 

(8) 
0) 



where j — 1,2,3, W is the Lorentz factor. The equa- 
tions 01 are closed by an equation of state (EOS) given 
by P = p(p 7 e). For an ideal gas, the EOS reads, 

p = (T - l)pe, (10) 

where T is the adiabatic index of the ideal gas. 

3. NUMERICAL SCHEMES 

To numerically solve Eq. each spatial dimension is 
discretized into cells. Using the method of lines, the time 
dependent evolution of Eq. 0] can be expressed in the 
semi-discrete form 



dt 



(11) 



Ax 

i,j + l/2,k r i,j-l/2,k 

Ay 

r i,j,k+l/2 F ij,k-l/2 



where Ff ±1/2j . fc , F» J . ±1/2 k and F^. fc±1/2 are the fluxes 
at the cell interface. 
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In order to achieve high-order accuracy in time, the 
time integration is done using a high-order t otal variation 
dimin ishing (TVD) Runge-Kutta scheme ijShu fc Osherl 
1988), which combines the first-order forward Euler steps 
and involves prediction and correction. For example, the 
third-order accuracy can be achieved via 

U«=U n + At£(U n ) (12) 
u( 2 > = ~U" + ~U« + ~AtL(U«) (13) 

U n+1 = ^U" + |u^ + |a£L(U< 2 )), (14) 
o o o 

where £(U) is the right hand side of Eq. El U™ +1 is 
the final value after advancing one time step from U™. 
Besides the third-order Runge-Kutta method (RK3), 
the standard fourth-order Runge-Ku tta (RK4) and fifth- 
order Runge-Kutta method (RK5) (Lambert 19911) can 
also be used for time integration. We usually use TVD 
RK3 in our calculations unless otherwise stated. We have 
also implemented and tested RK4 and RK5 and use them 
in some calculations. 

Note the method of lines, treats information from cor- 
ner zones differently than methods using time-averaged 
fluxes which make some use of corner information 
when implemented with Strang splitting. Recently, 
dimensionally unsplit methods have been developed 
l|Mignone. Plewa fc Bodoll20M ICardiner fc Stondl200l 
which use the corner-transport upwind (CTU) method of 
l)Colellalll990|) to include corner information completely. 

RAM, however, uses the third-order Runge-Kutta al- 
gorithm for time integration consisting of several forward 
Euler substeps which are corrected to achieve high accu- 
racy. Since the final results after one full time step are 
computed through several substeps, the corner cells (e.g., 
(xi+i, Uj+i, Zk+i)) of a zone at {xi,yj,Zk) affect its evo- 
lution in a full step, although this is not the case for an 
individual substep (§EH& §EH>- In this sense, RAM, 
is dimensionally unsplit. Comparison with codes using 
CTU methods on multi-dimensional test problems would 
be of value to perform in the future. We note that multi- 
dimensional test problems we have run with RAM pre- 
serve symmetries present in the initial conditions (see, 
e.g., Fig.iJl. 

To solve the fluxes across the cell interfaces, we have 
implemented two schemes. Both schemes involve a re- 
construction step to compute the variables at the cell 
interfaces. In the first class of schemes (§ I3.1fl . the re- 
construction is carried out on fluxes. The interface flux 
is obtained by, 

F m/2 = F(F i _ r! ...,F i+s ), (15) 

where the stencil (z — r, i — r + 1, + s) depends 
on the choice of the reconstruction scheme. This class 
can be considered finite difference schemes. We use 
F-X to denote this class of schemes, where X stands 
for the reconstruction scheme. In the other class of 
schemes (§ I3.2fl . the reconstruction is carried out on 
the unknown variables. Then the interface unknowns, 
U7^7, 2 = U (U;_ r , Uj +S ), are used to obtain interface 
fluxes, 

F m/2 = i?(U- 1/2 ,U+ 1/2 ), (16) 

where — and + denote the left and right side of the in- 
terface i + 1/2, respectively. The flux function R is an 



approximate Riemann solver. The cell unknowns Ui are 
considered as cell averages. Thus, the second class can 
be considered as finite volume schemes. We use U-X to 
denote the second class of schemes, where X stands for 
the reconstruction scheme. 

3.1. F-X Schemes: Reconstruction Of Fluxes 

In the F-X schemes, we have implemented the weighted 
essentially non-oscillatory (WENO) scheme and piece- 
wise linear method (PLM) for the reconstruction of 
fluxes. 

Essentia lly non-oscillator y (ENO) finite difference 
schemes fHarten et al.llT987Tl for equations of hyperbolic 
conservation laws use adaptive stencils in calculating the 
fluxes across the cell interfaces so that high-order ac- 
curacy can be achieved and numerical oscillations near 
discontinuities can be significantly reduced. Later to- 
tal variation diminishing (TVD) Runge-Kutta methods 
were applied to ENO schemes to make the schemes 
more computationally efficient, and the ENO reconstruc- 
tions were carried out on fluxes instead of cell averages 
l|Shu fc Osherlli 988. 1983). EN schemes have been fur- 
ther improved by many researchers (see iShul 119971 for 
a review). One improved scheme is the fourth-order 
weighted essentially non-os cillatory (WENO) schem e 
which was first introduced bv lLiu. Osher fc Chanl Jl994). 
In the WENO schemes, a weighted combination of sev- 
eral possible stencils are used instead of just one stencil. 
This would improve the accuracy while keeping the es- 
sentially non-oscillatory property near discontinuities. A 
new algorithm for computing the weights has been used 
in a modifie d fifth -order WENO scheme developed by 
l.liang fc ShulllT9"96T) . The WENO scheme oflJiang fc ShvJ 
l)1996l) has been applied bv lFeng et alJ l|2004l) to cosmo- 
logical simulations. In this paper, we have implemented 
the fift h-order modified WENO scheme of Jiang fc ShiJ 
( 1996) in special relativistic hydrodynamics. The WENO 
schem es for hyp erbolic conservation laws have various 
forms <|Shulll997|) . The particular WENO scheme we use 
is the characteristic-wise flux splitting finite difference 
scheme. 

In the WENO scheme, the partial differential equations 
arc discrctized into cells in spatial dimensions (Eq. 12). 
The evolution of the partial differential equations is 
solved by using a TVD Runge-Kutta scheme. The key 
component of the WENO scheme is to compute the fluxes 
at the cell interfaces. As an example, we shall consider 
the x-direction flux across the cell interface between xi 
and Xi+x. We shall assume that states at Xi-2, 
Xi+2, and Xi+3 are either in the computational domain 
or can be supplied by boundary conditions. 

First, we construct a Roe type average state at the cell 
interface x i+1 / 2 from the left state at xi and right state 
at Xi+i. The eigenvectors of the average state will be 
used as the basis for the decomposition of fluxes, F m , 
where m = i — 2, i — 1, i, i + 1, i + 2,i + 3. We use 
yfph as weight (|Eulderink fc Mellemalll995|) to compute 
the weighted averaged of pressure, density, and velocity. 
Strictly sp eaking, thi s does not satisfy all of the Roe's 
conditions l)Roelll98"H) . However, we do not use the av- 
erage state to compute the flux directly. Therefore strict 
Roe's conditions are not required. Indeed, simply taking 
arithmetic averages of the primitive variables also works 
very well. 
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For special relativistic hydrodynamics equations, the 
complete characteristic structure of these hy p erboli c 
equations have been derived by iDonat et al.l J.1998). 
Thus it is feasible to extend the characteristic-wise 
WENO to special relativistic hydrodynamics knowing 
the eigenvalues and eigenvectors of the Jacobian matri- 
ces of the equations. The eigenvalues of the Jacobian 
matrices, 

9F ^ (17) 



B 



<9U 



namely the speed of characteristic waves, can have dif- 
ferent signs. In other words, these waves could prop- 
agate towards different directions. The second step of 
the WENO scheme is to split the flux into left-going and 
right-going fluxes so that we can treat them separately. 
This flux splitting approach can avoid entropy violation 
and make the scheme more robust. We use the local 
Lax-Freiderichs splitting given by, 

F+=F + aU (18) 

F = F aU, (19) 

where a is the local maximum of the absolute values of 
wave speed for states at Xi-2,i-i,M+M+2,i+3- The eigen- 
values are all positive for the right-going flux F + , whereas 
they are negative for the left-going flux F . 

Because of the obvious symmetry between the left- 
going flux F - and right-going flux F + , we shall only con- 
sider here the right-going waves. The stencil for the right- 
going waves at i + 1/2, which consists of Xi—2,i— i,i,i+i,»+2 
in the fifth-order WENO scheme, is slightly upwind- 
biased. We can expand F+, where m = i — 2,i — 
1, i, i + l,i + 2, in terms of the five right eigenvectors 
ijDonat et al.l lT998) of the average state at i + 1/2, R„, 
where n = 1,2, 3, 4, 5. The expansion is given by 



^ ^ Cmn f^n j 



(20) 



where, m = i — 2, i— 1, i, i + 1, i + 2. It is very straightfor- 
ward to compute c rnn since we also know the left eigen- 
vectors, L„, of the Jacobian matrix. c rnn is given by, 

Cmn -^n ' F m . (21 ) 

Given c mn , we can construct the coefficients c i+1 / 2 .n 
for the cell interface at 2^+1/2 by giving different weights 
to different c mn . The weights depend upon the smooth- 
ness of the stencils, and smooth stencils are given more 
weight. This results in the fifth-order accuracy in smooth 
region and ENO near discontinuities. For the details of 
how the weights are chosen in the WENO r econstruction 
schem e, we refer the readers to the work of Uiang fc Shul 

Using the coefficients Cj+i/2, n ) where n denotes char- 
acteristic waves, the right-going flux at the cell interface 
at x i+1/2 is given by 



F+ 

1+1/2 



c i+l/2,nR-n- 



(22) 



Because of the obvious symmetry, the left-going flux, 
1^1+1/2' can be computed in the similar procedure. Note 
that the left-going flux is also upwind-biased, namely 
right-biased. 



In computing the WENO weights, a parameter, e, is 
introduced to avoid denominator being zero. This is the 
only free parameter in the WENO scheme. Moreover, 
the results are insensitive to the value of e as long as it 
is a small number about 10 -6 . 

We use F-WENO to denote the above method, which 
uses the WENO scheme for the characteristic-wise re- 
construction of splitted fluxes. The reconstruction of the 
coefficients Cj+i/2, n from c mn can also be computed using 
other high-order reconstruction alg orithms instead of th e 
WENO reconstruction algorithm of Jian g fc Shul l)1996j) . 
Besides the WENO, we have also implemented the PLM 
for the reconstruction with a gener alized minmod slope 
limiter ijKurganov fe Tadmorll2000j) . Given Ci-\, Cj, and 
the left-biased interface value reads, 

c i+i/2 = c i + 0.5minmod(#(ci — Ci-x), (23) 

0.5(cj+i - Ci-i),0{ci+x - Ci)), 

where 1 < 6 < 2, and the minmod function reads, 

minrnod(x, 2/, z) = — (sgn(x) + sgn(y)) (24) 
(sgn(x) +sgn(z))min(|x|, \y\, \z\), 

here the sgn function returns the sign of the number. 
This becomes the more diffusive normal minmod lim- 
iter when 6 = 1, an d becomes the m onotonized central- 
difference limiter of Ivan Leerl llT977h when 6 = 2. We 
usually use 6 = 1.5 unless otherwise stated. We call 
this method F-PLM. Our numerical tests show that the 
results of F-PLM are comparable to those of F-WENO 
(§0- 

3.2. U-X Schemes: Reconstruction Of Unknowns 

Besides the F-X schemes (§ I3.1fl . we have also im- 
plemented another method of computing the flux at 
the interface to work with the Runge-Kutta scheme 
for time integration of Eq. El In this class of 
schemes, U-X, instead of reconstructing flux directly 
as in the F-X schemes, a left state and a right state 
at the interface are reconstructed, and then the flux 
at the interface is calculated using these two states. 
A lot of recent methods, including the so-called high- 
resolution shock-ca pturing (HRSC) methods like the 
GENESIS method llAlov et all I1999D and some high- 
resolution central schemes {Lucas- Serrano et al.l 120041: 
iDel Zanna &; Bucciantinill2002|) . can be considered to be 
in this category. Their main difference is how to compute 
the flux at the interface given the left and right states 
and how to reconstruct interface states. In the HRSC 
methods, an approximate Riemann solver utilizing the 
information of the characteristic waves is used, whereas 
the central schemes use simplified expressions for the flux 
without using the characteristic information beyond the 
fastest wave speed, which is nevertheless required for the 
Courant-Friedrich-Levy (CFL) condition. 

The reconstruction of states at the cell interfaces 
is usually performed by a high-order interpolation 
method like the piecewise parabolic method (PPM) 
llColella fc Woodwardl fl984) or piecewise linear method 
(PLM) . We have implemented both reconstruction meth- 
ods in the RAM code. For the parameters in the 
PPM algorithm, we choose the values suggested by 
iMarti fc Mulled <jl996j) for almost all the tests with PPM. 
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We consider it undesirable to fine-tune these parameters 
to achieve better agreement with analytic solutions for 
specific tests. Doing so, we believe, can create a false 
sense of how a given code performs in general simulations 
for which the exact solution is not previously known. 
For the PLM algorithm, we use a generalized minmod 
slope limiter as described in § 13.11 Again, the parameter 
6 = 1.5 is used by default. The reconstruction is per- 
formed on the so-called primitive variables instead of con- 
served variables because unphysical conserved variables 
may arise otherwise. The pressure and the proper den- 
sity are reconstructed directly, whereas velocities are re- 
constructed using a combination of reconstructing three- 
velocity and Lorentz factor. Since both the velocities and 
Lorentz factor are reconstructed, they are unlikely to be 
self consistent. To make them consistent, four-velocities 
are derived by multiplying three-velocities with Lorentz 
factor. Using these four-velocities, new consistent veloc- 
ities and Lorentz factor could be derived. We found that 
this procedure is usually more robust than using only 
three-velocities or only four-velocities for the construc- 
tion of velocity and Lorentz factor. 

To compute the flux across the cell interface given the 
reconstructed right and left states at the interface, we 
have implemented several Riemann solvers i ncluding the 
modified Marquina's flux ijAlov et al] 11999]) . local Lax- 
Fried richs flux and relativistic HLLE IjSchneide r et alJ 
1993) in the RAM code. A comparison of several 
schemes for computing the flux has been performed by 
ILucas-Serrano et al.l ((2004) . They have shown that the 
results are relatively independent of the flux schemes. 
For the results shown in this paper, the modified Mar- 
quina's algorithm is used to compute the interface fluxes. 

3.3. Failsafe Time Integration and Root Finder 

Since our code is explicit, the time step for integration 
is subject to the Courant-Friedrich-Levy (CFL) condi- 
tion. We usually choose a CFL number to be 0.2-0.5. 
Unfortunately, it is not always failsafe. Occasionally 
unphysical results can be produced for ultra-relativistic 
flows unless a very small CFL number or a diffusive 
scheme is used. But a very small CFL number is some- 
times unacceptable because of the time it costs. To make 
the code robust, a fall-back mechanism is employed in 
our code, though we found that such failures rarely oc- 
cur. That is a normal CFL number is used for the first 
try of time integrations. If it fails, the code will return 
to the beginning of the failed step and use more diffu- 
sive schemes for reconstruction and/or use a smaller time 
step. The recalculation with more diffusive schemes is 
performed on the whole grid, when AMR is not being 
used. The fall-back method on AMR grids will be fur- 
ther discussed in § [3] The fall-back requires us to keep 
a copy of the original states of conserved variables for 
a period of the whole Runge-Kutta step in order to be 
able to fall back. Fortunat ely, the Runge-Kutta scheme 
we use llShu fc Osherll98 8) requires the original states to 
be saved anyway (Eqns.ElEU Thus it does not 

increase the memory usage of the computation. In the 
subsequent steps, the normal reconstruction scheme will 
be used and the CFL number will gradually increase to 
its normal value if it was previously decreased. The fall- 
back mechanism makes the code almost failsafe. If un- 
physical results persist even when the first-order method 



is employed and a very small time step is used, the nu- 
merical simulation will stop. One possible remedy when 
the code runs out of methods is to set the pressure of 
bad cells to a floor value. Fortunately this has never 
happened in our simulations. 

In each time step, conserved variables U = 
(D, S 1 , S 2 , S 3 , t) t are updated directly. In order to cal- 
culate the fluxes at the cell interfaces, physical variables: 
pressure, proper rest mass density and velocity need to 
be computed. The processes of recovering physical vari- 
ables from conser ved variables involve a qu artic equation 
for velocity (e.s.. iDuncan fc Hug-hesllT994T) . Though the 
quartic can be solved analytically, computing the ana- 
lytic solution is actually more expensive than a simple 
numerical root finder, such as Newton-Raphson itera- 
tion. Before the Newton-Raphson iteration, a certain 
condition, (r + D) 2 > D 2 + S 2 , is checked to make sure 
that physical results can be produced from given con- 
served variables. Sometimes, unphysical results, such as 
negative pressure and large velocity v > 1, will be pro- 
duced. Then the fall-back mechanism will be used until 
physical results are produced or the maximum allowed 
number of fall-back is reached. 

4. TEST PROBLEMS 

For any numerical code, it is very important to do 
substantial tests. We have done a series of tests with 
our code. Some of the tests are the so-called Riemann 
problems, which consist of computing the decay of an 
initial discontinuity of two constant states. It is possible 
to get exact solut ions for relativistic Riemann problems 
ijPons et al.ll2000() . We have also done some extensively 
studied problems which have no analytic solutions. In 
all cases, our results are comparable to published results. 
In this section, we will present our numerical results for 
some standard tests. We will compare four schemes (§E}: 
the finite difference characteristic- wise WENO, the fi- 
nite difference characteristic-wise PLM, the finite volume 
component-wise PPM and the finite volume component- 
wise PLM, denoted by F-WENO, F-PLM, U-PPM, and 
U-PLM, respectively. A CFL number of 0.5 is used for 
these tests unless otherwise stated. 

4.1. One- Dimensional Riemann Problem 1 

In this test, the one-dimensional numerical region 
(0 < x < 1) initially consists of two constant states: 
p L = 13.33, p L = 10.0, v L = 0.0 and p R = 10~ 8 , 
Pr = 1.0, vr = 0.0, where L stands for the left state, and 
R the right state. The fluid is assumed to be an ideal 
gas with an adiabatic index T — 5/3. The initial discon- 
tinuity is at x = 0.5. The results are shown in Fig. ^ 
In this test problem, the evolution of the initial discon- 
tinuity gives rise to a shock, a rarefaction wave, and a 
contact discontinuity. This is a fairly easy test. All mod- 
ern special relativistic hydrodynamics codes can capture 
the expected features, acquire correct positions of the 
shock front, contact discontinuity and rarefaction wave. 
To quantitatively measure the errors, we calculated the 
L\ norm errors, L\ = Axj\uj — u(xj)\, where u(xj) 
is the exact solution at Xj and Uj is the numerical re- 
sult. The L\ norm errors of density for four schemes 
with various grid resolutions at t = 0.4 are shown in Ta- 
ble[T] The accuracy of our results is comparable to t hat of 
ILucas-Serrano et all (|2004l) : lMarti fc Mullerl (fl996). The 
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Fig. 1. — One-dimensional Riemann problem 1 at t = 0.4. 
Results for four schemes: (a) F-WENO; (b) F-PLM; (c) U-PPM; 
(d) U-PLM are shown. The computational grid consists of 400 
zones. Numerical results are shown in symbols, whereas the exact 
solution is shown in solid lines. We show proper mass density 
(square), pressure (triangle) and velocity (plus sign). 



Fig. 2. — One-dimensional Riemann problem 2 at t = 0.4. 
Results for four schemes: (a) F-WENO; (b) F-PLM; (c) U-PPM; 
(d) U-PLM are shown. The computational grid consists of 400 
zones. Numerical results are shown in symbols, whereas the exact 
solution is shown in solid lines. We show proper mass density 
(square), pressure (triangle) and velocity (plus sign). 



TABLE 1 

li errors of the density for the id rlemann 

Problem 1. Four schemes with various 
resolutions with uniform spacing are shown 

AT t = 0.4. 



Scheme 


7V a 


Li Error 


Convergence Rate 


F-WENO 


100 


1.31e-l 






200 


7.25e-2 


0.85 




400 


3.32e-2 


1.1 




800 


2.08e-2 


0.67 




1600 


1.00e-2 


1.1 




3200 


5.07e-3 


0.98 


F-PLM 


100 


1.47e-l 






200 


8.50e-2 


0.79 




400 


4.06e-2 


1.1 




800 


2.33e-2 


0.80 




1600 


1.22c-2 


0.93 




3200 


7.48e-3 


0.71 


U-PPM 


100 


1.27c-l 






200 


7.30e-2 


0.80 




400 


3.47e-2 


1.1 




800 


1.97e-2 


0.82 




1600 


9.77c-3 


1.0 




3200 


5.10c-3 


0.94 


U-PLM 


100 


1.32e-l 






200 


8.57e-2 


0.62 




400 


3.86e-2 


1.2 




800 


2.27e-2 


0.77 




1600 


1.15e-2 


0.98 




3200 


6.48e-3 


0.83 



a Number of grid points 



order of the convergence rate is about 1. This is consis- 
tent with the fact that there are discontinuities in the 
solution. 

4.2. One-Dimensional Riemann Problem 2 

In this test, the one-dimensional numerical region 
(0 < x < 1) initially consists of two constant states: 
p L = 1000.0, p L = 1.0, v L = 0.0 and p R = 10~ 2 , 
p R = 1.0, v R = 0.0, where L stands for the left state, 
and R the right state. The fluid is assumed to be an 
ideal gas with an adiabatic index T — 5/3. The initial 
discontinuity is at x = 0.5. The results are shown in 



Fig. [3 In this test problem, the evolution of the ini- 
tial discontinuity gives rise to a right-moving shock, a 
left-moving rarefaction wave, and a contact discontinu- 
ity in between. Behind the shock, there is an extremely 
thin dense shell. The width of the shell is only 0.01056 
at t — 0.4. With 400 uniform zones for < x < 1, 
the thin shell is only covered by 4.2 zones. Due to the 
smearing at the contact discontinuity and the shock, it 
is not surprising that the thin shell is not well resolved 
in our results with only 400 zones. At this resolution 
the maximum density in the shell for F-WENO, F-PLM, 
U-PPM and U-PLM is 79%, 69%, 72% and 63% of the 
analytic value. However, there is no difficulty in resolv- 
ing the thin shell with increased resolution. We have 
calculated the L\ norm errors of density for various nu- 
merical resolutions (see Table 0) and achieve convergence 
as expected for problems with sharp discontinuities. The 
results are consistent with other published results (e.g., 
iLucas-Serrano et al1l2004f) . 

4.3. One-Dimensional Riemann Problem 3 

In this test, the one-dimensional numerical region (0 < 
x < 1) initially consists of two constant states: pl = 1-0, 
PL = 1.0, v L = 0.9 and p R = 10.0, p R = 1.0, v R = 0.0, 
where L stands for the left state, and R the right state. 
The fluid is assumed to be an ideal gas with an adiabatic 
index T = 4/3. The initial discontinuity is at x = 0.5. 
The results are shown in Fig. [21 In this problem a strong 
reverse shock forms in which post-shock oscillations are 
visible for the U-PPM and U-PLM simulations, espe- 
cially in the pressure profile. In the F-WENO simulation, 
the post-shock pressure oscillations is smaller those in 
the U-X simulations. In the F-PLM simulation the post- 
shock pressure oscillation is nearly invisible to the eye, 
though they are still present at the 0.1% level. Reducing 
the CFL number from 0.5 decreases the post-shock oscil- 
lations but does not eliminate them completely Table |31 
presents the L\ norm errors and L\ order convergence 
rates for this problem. Note that the U-PPM scheme 
is more accurate and converges faster for this problem 
than other schemes. This is due to the fact that U-PPM 
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TABLE 2 

l\ errors of the density for the id rlemann 

Problem 2. Four schemes with various 
resolutions with uniform spacing are shown 

AT t = 0.4. 



Scheme 


N 


L\ Error 


Convergence Rate 


F-WENO 


100 


2.10e-l 






200 


1.42e-l 


0.56 




400 


9.29e-2 


0.61 




800 


5.54e-2 


0.75 




1600 


2.54e-2 


1.1 




3200 


1.51e-2 


0.75 


F-PLM 


100 


1.96c-l 






200 


1.42e-l 


0.46 




400 


1.06e-l 


0.42 




800 


7.21e-2 


0.56 




1600 


3.92e-2 


0.88 




3200 


2.44e-2 


0.68 


U-PPM 


100 


2.18e-l 






200 


1.52e-l 


0.52 




400 


9.52e-2 


0.68 




800 


5.42e-2 


0.81 




1600 


2.67e-2 


1.0 




3200 


1.67e-2 


0.68 


U-PLM 


100 


2.13c-l 






200 


1.65e-l 


0.37 




400 


1.25e-l 


0.41 




800 


8.68e-2 


0.53 




1600 


4.49e-2 


0.95 




3200 


2.71e-2 


0.73 



a Number of grid points 



TABLE 3 

l\ errors of the density for the id rlemann 
Problem 3. Four schemes with various 
resolutions are shown at t = 0.4. 



OCilGmG 


N a 




vonvcr^i m k (' \. \ cit't 1 


F-WENO 


100 


9.97e-2 






200 


6.29e-2 


0.67 




400 


3.01e-2 


1.1 




800 


1.69e-2 


0.83 




1600 


9.48e-3 


0.83 




3200 


5.24e-3 


0.86 


F-PLM 


100 


1.12c-l 






200 


6.98e-2 


0.68 




400 


3.45e-2 


1.0 




800 


1.94e-2 


0.83 




1600 


1.13c-2 


0.78 




3200 


6.54e-3 


0.79 


U-PPM 


100 


9.72c-2 






200 


5.60e-2 


0.80 




400 


2.49e-2 


1.2 




800 


1.30e-2 


0.94 




1600 


6.06e-3 


1.1 




3200 


3.11c-3 


0.96 


U-PLM 


100 


9.53c-2 






200 


6.32e-2 


0.59 




400 


2.99e-2 


1.1 




800 


1.78e-2 


0.75 




1600 


1.04e-2 


0.78 




3200 


6.10e-3 


0.77 



a Number of grid points 
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Fig. 3. — One-dimensional Riemann problem 3 at t = 0.4. Re- 
sults for four schemes: (a) F-WENO; (b) F-PLM; (c) U-PPM; (d) 
U-PLM arc shown. The computational grid consists of 400 zones. 
Numerical results are shown in symbols, whereas the exact solu- 
tion is shown in solid lines. We show proper mass density (square), 
pressure (triangle) and velocity (plus sign). This problem contains 
a strong reverse shock in which post-shock pressure oscillations are 
visible in panels (c) and (d). 



scheme captures the contact discontinuity in fewer zones 
©■ 

4.4. One-Dimensional Riemann Problem With 
Non-Zero Transverse Velocity: Easy Test 

Many problems of interest in hydrodynamics involve 
strong shear flows. Astrophysical jets include shearing 
layers where mixing of ambient material into the fast jet 
flow may be important. It is therefore important to test 
the ability of numerical codes to handle Riemann prob- 



lems with velocity components transverse to the direction 
of propagation of the main flow. In Newtonian hydro- 
dynamics the transverse momentum is simply advected 
with the flow and is not coupled directly to the equa- 
tion of motion in the longitudinal direction. No serious 
difficulty is introduced by the presence of transverse ve- 
locity components though transverse kinetic energy dis- 
sipated through viscous dissipation can alter the longitu- 
dinal motion by affecting the pressure. However, in rela- 
tivistic flow, transverse velocities are directly coupled to 
the dynamics along all directions by the Lorentz factor 
which depends on all velocity components. This cou- 
pling makes relativistic Riemann problems with trans- 
verse velocity much more difficult to solve correctly in 
a numerical code. In particular, higher resolution is 
needed. Under-resolved simulations produce incorrect 
shock positions along the normal direction. Here, we 
present relatively easy ID tests of Riemann problems in- 
cluding transverse velocity which can be resolved with 
moderate resolution. In § I6.ll we show tests requiring 
high resolution and demonstrate the ability of adaptive 
mesh refinement to accurately simulate very relativistic 
flows with significant transverse velocity components. 

In this test, the one-dimensional numerical region 
(0 < x < 1) initially consists of two constant states: 
PL = 1000.0, PL = 1.0, v xL = 0.0, v yL = 0.0 and 
PR = 10~ 2 , pr = 1.0, v xR = 0.0, v yR = 0.99, where 
L stands for the left state, and R the right state. The 
fluid is assumed to be an ideal gas with an adiabatic in- 
dex r = 5/3. The initial discontinuity is at x = 0.5. 
The results are shown in Fig. 0] In this test, the pres- 
ence of transverse velocity alters the structure of the 
Riem ann problem ijPons et alJ l2000: Rcz zolla &: Zanottl 
2002). This problem is relatively easy because the trans- 
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Fig. 4. — "Easy" one-dimensional Riemann problem with non- 
zero transverse velocity (v y = 0.99c) at t = 0.4. Results for four 
schemes: (a) F-WENO; (b) F-PLM; (c) U-PPM; (d) U-PLM are 
shown. The computational grid consists of 400 zones. Numerical 
results are shown in symbols, whereas the exact solution is shown 
in solid lines. We show proper mass density (square), pressure 
(triangle) and velocity in x-direction (plus sign). 
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Fig. 5. — One-dimensional shock heating problem in planar 
geometry at t = 2.0. Results for four schemes: (a) F-WENO; (b) 
F-PLM; (c) U-PPM; (d) U-PLM are shown. The parameter 8 in 
the minmod slope limiter for U-PLM is set to be 1.0 in this test. 
The computational grid consists of 100 zones. Numerical results 
are shown in symbols, whereas the exact solution is shown in solid 
lines. We show proper mass density (square), pressure (triangle) 
and velocity (plus sign). 



TABLE 4 

li errors of the density for the "easy" id 
Riemann problem with non-zero transverse 

velocity at t = 0.4. four schemes with 
various resolutions using a uniform grid are 

SHOWN. 



Scheme 


N a 


Li Error 


Convergence Rate 


F-WENO 


100 


7.58e-l 






200 


3.92e-l 


0.95 




400 


2.31e-l 


0.76 




800 


1.18c-l 


0.97 




1600 


6.58e-2 


0.84 




3200 


3.44e-2 


0.94 


F-PLM 


100 


8.26c-l 






200 


4.59e-l 


0.85 




400 


2.77c-l 


0.73 




800 


1.49e-l 


0.89 




1600 


8.00e-2 


0.90 




3200 


4.63e-2 


0.79 


U-PPM 


100 


8.48e-l 






200 


4.25e-l 


1.0 




400 


2.41e-l 


0.82 




800 


1.27e-l 


0.92 




1600 


6.43e-2 


0.99 




3200 


3.34e-2 


0.95 


U-PLM 


100 


9.00c-l 






200 


4.72e-l 


0.93 




400 


2.88e-l 


0.71 




800 


1.52e-l 


0.92 




1600 


8.86e-2 


0.78 




3200 


4.95e-2 


0.84 



a Number of grid points 



verse velocity is in the cold gas of the right state, not in 
the relativistically hot left state or in the rarefaction fan 
which subsequently propagates into it. The simulation 
agrees well with the analytic solution even at the rela- 
tively modest resolution of 400 zones. Table 01 presents 
the L\ norm errors and L\ order convergence rates for 
this problem. 

4.5. One-Dimensional Shock Heating Problem In 
Planar Geometry 



The shock heating problem is a standard test to study 
the ability of a code to handle very strong shocks with 
sufficiently few zones and without excessive post shock 
oscillations. Cold gas flows toward a reflecting boundary 
at x — 1.0 and a reverse strong shock forms and propa- 
gates to the left decelerating the gas to zero speed. The 
gas has an initial speed of v = 0.9999999999 = 1.0-10" 10 
(corresponding to a Lorentz factor of W — 70710.675) 
and initial density of p = 1.0. Due to conservation of en- 
ergy and the fact the gas initially has nearly zero internal 
energy compared to kinetic energy (we use a small value 
eo = 0.003 for numerical reasons), the specific internal 
energy after the shock is simply e = W — 1. The com- 
pression ratio across the shock can be arbitrarily large in 
the relativistic case and is given by 



where V is the adiabatic index of the gas. In this case 
T = 4/3 so (7 = 282845.70 and the shock velocity is given 

by v s = {T ~w ) +i H = 0.33332862. In Fig. we show our 
results for a uniform mesh of 100 zones compared with 
the analytic solution at t = 2.0. 

We note that for this problem with a constant state 
behind the shock, as opposed to a thin shell which might 
more naturally occur in a realistic flow, the maximum 
Lorentz factor that can be achieved numerically is just 
limited by floating point precision. 

The parameter 9 in the minmod slope limiter for U- 
PLM is set to be 1.0 in this test to eliminate strong oscil- 
lations which would appear if the default value 8 = 1.5 is 
used in U-PLM. In this test with U-PPM, stronger dissi- 
pation is also needed to avoid crashes. In particular, one 
of the PPM parameters, e^ 2 \ is s et to be 10~ 4 instead 
of the default value = 1.0 fsee iColella fc Woodwardl 
IT9841 for the meaning of this parameter). The reflecting 
wall at x — 1.0 poses difficulties for numerical simula- 
tions and gives rise to visible errors for zones near the 
wall. The errors of density at the nearest zone to the 
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reflecting wall are 3.9%, 2.4%, 7.4%, and 4.3%, for F- 
WENO, F-PLM, U-PPM, and U-PLM, respectively. For 
this particular problem with very strong shock but sim- 
ple structure, more diffusive schemes F-PLM and U-PLM 
perform better than F-WENO and U-PPM. 

4.6. Isentropic Smooth Flows 

Most tests presented in this paper involve strong 
shocks. In addition to tests with discontinuities, it is 
also very important to test the capability of the schemes 
to handle smooth flows. In this section, we present two 
tests involving the isentropic evolution of smooth flows 
in ID and 2D. 

In the first test, a one-dimensional isentropic smooth 
pulse is set up in a uniform reference state, and the run 
stops some time before a shock is formed. T he test is sim- 
ilar to the convergence tests performed bv iColella et al J 
(2006) for their Newtonian hydrodynamics code. The 
initial density structure at t = is given by 

p (x) = Pici (l + af(x)), (26) 

where p lc f is the density of the reference state and 



/w = U WL)2 " 1)4 



H < L (27) 
otherwise, v ' 



here a is the amplitude of the pulse and L the width 
of the pulse. The pressure is given by the isentropic 
relation, p — Kp r , where K is a constant. The initial 
velocity of the reference state is v TC { = 0. The initial 
velocity inside the pulse at t = is set up by assuming 
one of the two Riemann invariants, 



-ln(- 
2 1 



1 



ln( 



(28) 



is constant across the whole region, where c s is the sound 
speed. We note that the other Riemann invariant, 



1 



ln( 



J+ = -ln( ) -I 



(29) 



is not constant. The exact solution of the test can be 
obtained by using standard characteristic analysis. The 
pulse will have a smooth shape before a shock eventu- 
ally forms. The width and height of the pulse does not 
change before the shock forms. But the pulse will be- 
come increasingly asymmetric as the shape of the front 
of the pulse becomes steeper during the propagation (see 
Fig. EJ) . Behind the moving pulse, the fluid goes back to 
the reference state. 

Our computational region for this test is —0.35 < x < 
1. The reference state is, j? re f = 100, p le { = l,i> re f = 0, 
the amplitude of the pulse is a = 1.0, and the width is 
L = 0.3. The adiabatic index in the equation of state is 
r = 5/3. 

We have run this test with different schemes and var- 
ious numerical resolution. Because the flow is very 
smooth, all methods perform very well for this test. 
Fig. El shows the exact solution and the numerical re- 
sults of F-WENO with 80 uniform zones. The results 
of the convergence rates are shown in Table |3J Various 
Runge-Kutta methods, including the third-order TVD 
scheme, standard fourth-order scheme (RK4) and fifth- 
order scheme (RK5) have been used in this test. We find 




Fig. 6. — One-dimensional isentropic flow. The initial pulse is 
shown on the left, and the pulse at t = 0.8 on the right. Pressure, 
density and velocity are shown in the top, middle and bottom 
panels respectively. The solid lines are exact solutions. Numerical 
results of F-WENO at t = (plus signs) and t = 0.8 (triangle) are 
shown. 



that the F-WENO scheme, which is fifth-order accurate 
in space and third-order accurate in time, has an L\ or- 
der convergence rate of ~ 3 — 4 for this one-dimensional 
test with smooth flow. For the F-WENO-RK4 and F- 
WENO-RK5 schemes, which are fourth and fifth-order 
accurate in time integration respectively, order of con- 
vergence rates up to ~ 5 can be achieved. Note that 
the difference between the results of RK4 and RK5 is 
very small. For schemes other than WENO, including 
schemes using RK4 or RK5, the order of convergence 
rate is about 2 (Table |5j). These results clearly show 
that the fifth-order WENO scheme is very accurate and 
converges very quickly for smooth flows. 

PPM is generally 3rd-order accurate in space, but the 
flattening procedure and the requirement of monotonic 
profiles degrades the accuracy of the scheme at places like 
local extrema or where the second derivatives of variables 
are large. Thus its accuracy can be as low as first-order 
in some places. Moreover, in the U-PPM scheme, the 
reconstruction is carried out on primitive variables, not 
conservative variables. The conservative variables are 
cell averages. But the primitive variables converted from 
conservative variables are not exactly cell averages. How- 
ever, the PPM reconstruction algorithm regards them as 
cell averages. Therefore the interface values of primitive 
variables may not have third-order accuracy. Thus it is 
reasonable that the U-PPM scheme does not achieve a 
third-order convergence rate. 

In the second test, we have performed two-dimensional 
calculations to assess the convergence rate of the WENO 
scheme in multi-dimension. The computational region in 
this test consists of a two-dimensional box in Cartesian 
coordinates with 0.0 < x < 3.75 and 0.0 < y < 5.0. The 
boundary conditions are periodic for all four sides of the 
box. Like the one-dimensional test of isentropic flows, 
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TABLE 5 

L\ ERRORS OF THE DENSITY FOR THE ID ISENTROPIC 
FLOW PROBLEM AT i = 0.8. RESULTS WITH VARIOUS 
RESOLUTIONS USING A UNIFORM GRID ARE SHOWN. 



TABLE 6 

L\ ERRORS OF THE DENSITY FOR THE 2D ISENTROPIC FLOW 
PROBLEM AT t = 2.4. RESULTS WITH VARIOUS RESOLUTIONS 
USING A UNIFORM GRID ARE SHOWN. 



Scheme 3, 


N° 


Li Error Convergence Rate 


F-WENO 


80 


2.07e-3 






160 


1.10e-4 


4.2 




320 


1.70e-5 


2.7 




640 


1.47e-6 


3.5 




1280 


1.58e-7 


3.2 




2560 


1.91e-8 


3.1 




5120 


2.37e-9 


3.0 


F-WENO-RK4 


80 


1.87e-3 






160 


1.18e-4 


1.0 




320 


1.31e-5 


3.2 




640 


6.80e-7 


4.3 




1280 


2.54e-8 


4.7 




2560 


7.91e-10 


5.0 




5120 


2.38e-ll 


5.1 


F-WENO-RK5 


80 


1.87e-3 






160 


1.17e-4 


1.0 




320 


1.30e-5 


3.2 




640 


6.82e-7 


4.3 




1280 


2.54e-8 


4.7 




2560 


8.01e-10 


5.0 




5120 


2.40e-ll 


5.1 


F-PLM 


80 


8.79c-3 






160 


4 05e-3 


1.1 




320 


1.22e-3 


1 7 




640 


3.10e-4 


2.0 




1280 


7.83e-5 


2.0 




2560 


1.96e-5 


2.0 




5120 


4.92e-6 


2.0 


F-PLM-RK4 


80 


8 85e-3 






160 


4 06e-3 


1.1 




320 


1.21e-3 


1.7 




640 


3.11e-4 


2.0 




1280 


7.84e-5 


2.0 




2560 


1.97e-5 


2.0 




5120 


4.93e-6 


2.0 


U-PPM 


80 


l.lle-2 






160 


2.47e-3 


2.2 




320 


7.02e-4 


1.8 




640 


1.38e-4 


2.3 




1280 


2.92e-5 


2.2 




2560 


6.48e-6 


2.2 




5120 


1.52e-6 


2.1 


U-PPM-RK4 


80 


1.10e-2 






160 


2.56e-3 


2.1 




320 


5.74e-4 


2.2 




640 


1.34e-4 


2.1 




1280 


3.10e-5 


2.1 




2560 


7.33e-6 


2.1 




5120 


1.82e-6 


2.1 


U-PLM 


80 


1.12e-2 






160 


3.56e-3 


1.7 




320 


1.03e-3 


1.8 




640 


2.61e-4 


2.0 




1280 


6.50e-5 


2.0 




2560 


1.62e-5 


2.0 




5120 


4.03e-6 


2.0 


U-PLM-RK4 


80 


1.12e-2 






160 


3.56e-3 


1.7 




320 


1.02e-3 


1.8 




640 


2.60e-4 


2.0 




1280 


6.49e-5 


2.0 




2560 


1.62e-5 


2.0 




5120 


4.04e-6 


2.0 


a RK4 and RK5 


denote the fourth 


and fifth-order 


Runge-Kutta methods, 


respectively. 


The third-order 



Runge-Kutta method (RK3) is used unless otherwise 
stated. 
b Number of grid points 



Scheme' 


A 


Li Error 


Convergence Rate 


F-WENO 


48 x 64 


7.35e-2 






96 x 128 


4.43e-3 


4.1 




192 x 256 


8.04e-4 


2.5 




384 x 512 


9.62e-5 


3.1 




768 x 1024 


1.12e-5 


3.1 


F-WENO- RK4 


48 x 64 


7.24e-2 






96 x 128 


4.75e-3 


3.9 




192 x 256 


4.70e-4 


3.3 




384 x 512 


3.18e-5 


3.9 




768 x 1024 


1.24e-6 


4.7 


F-WENO-RK5 


48 x 64 


7.19c-2 






96 x 128 


4.67e-3 


3.9 




192 x 256 


4.61e-4 


3.3 




384 x 512 


3.13e-5 


3.9 




768 x 1024 


1.22e-6 


4.7 



a The third, fourth (RK4), and fifth-order (RK5) Runge-Kutta 
methods are used with the fifth-order WENO scheme. 
b Number of grid points in x and y-direction 



there is a static uniform reference state, which is set to 
p I0 f = 100,/Jrcf = l,i>ref = 0. Pulses which are periodic 
in space are added to the system. Along the direction 
of k = (4/5,3/5), the profile is periodic with a spatial 
period of S — 3.0, and the profile is constant along the 
direction perpendicular to the vector k. Thus the spatial 
periods along the x and y-direction are 3.75 and 5.0, 
respectively. Note that these are consistent with the size 
of the computational box with periodic boundaries. The 
pulses move along the direction of the vector k. The 
initial density profile is given by po(d) (Eas. Eol fc ETjl . 
where d, the distance to the center of the nearest pulse, 
is given by d = mod(k • r + S/2, S) — S/2, here mod(a, b) 
returns the reminder of the division a/b, and r = (x,y). 
The amplitude of the pulse is a — 1.0, and the width is 
L = 0.9. The adiabatic index in the equation of state 
is T = 5/3. Similar to the one-dimensional test, the 
initial pressure is given by the isentropic relation, and 
the initial velocity by assuming the Riemann invariant, 
J_ is constant. 

We have run the two-dimensional test using the WENO 
scheme with three Runge-Kutta methods: RK3, RK4 
and RK5. Results of the L\ norm errors and convergence 
rate are shown in Tabled Similar to the one-dimensional 
test, both F-WENO-RK4 and F-WENO-RK5 performs 
better than F-WENO, which uses RK3. But the dif- 
ference of errors between F-WENO-RK4 and F-WENO- 
RK5 is very small. Therefore, for this test it is not worth 
using the more expensive RK5 for time integration. 

4.7. Two-Dimensional Tests: Wind Tunnel With Step 

In order to test the ability of our code to handle strong 
shocks in multiple dimensions, we have performed stan- 
dard tests published in the literature. 

Th e Emery step l|Emervl Il968t iWoodward fc Colellal 
Il984|) consists of a horizontal wind flowing into a step, 
represented as a reflecting boundary condition. The cor- 
ner of the step represents a singular point of the rar- 
efaction fan. As the wind collides with the step a reverse 
shock propagates back into the wind forming a bow shock 
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Fig. 7. — Emery step at t = 4.0 with 240 X 80 resolution. 
Thirty equally spaced contours of the logarithm of proper density 
are plotted. Results for four schemes: (a) F-WENO; (b) F-PLM; 
(c) U-PPM; (d) U-PLM are shown. 



which reflects off the upper boundary and forms a Mach 
stem which should remain straight and initially almost 
vertical. Since this structure can remain aligned with the 
vertical coordinate of the numerical grid, problems such 
as "odd-even decoupling" in some n umerical sche mes can 
act to incorrectly modify the stem <)Quirldll994l) . 

We simulate the Emery step in a computational box 
with < x < 3 and < y < 1 with 240 zones in the 
x-direction and 80 in the y-direction. The step is repre- 
sented as a reflecting boundary of height 0.2 beginning 
at x = 0.6 and continuing along the length of the box. 
The upper boundary and lower boundary for x < 0.6 
are also both reflecting. Initially the box is filled with 
p = 1.4 gas moving at v x = 0.999. The left boundary is 
an inflow with the same quantities. The right boundary 
is outflow. The adiabatic index in the equation of state 
is T = 1.4. The Newtonian Mach number is 3.0, and the 
corresponding relativistic Mach number M. — WM/W S , 
where M is the classical Mach number, W is the Lorentz 
factor of the gas and W s is the Lorentz factor of sound 
speed in the gas, is ~ 63. 

Our results are sh own in Fig. [71 These results are com- 
parable to those of iLucas-Serrano et al.1 l)2004() . Some 
numerical methods require special entropy fixes to accu- 
rately simulate the flow near the step corner. No such 
code modification was implemented or required for the 
simulations presented here because we are only interested 
in the global structures, especially the shocks. Though 
numerical boundary artifacts exist near the corner, the 
global solutions are not affected. 

The F-PLM and F-WENO methods are clearly supe- 
rior for this problem. U-PPM and U-PLM are noisy both 
downstream of the reverse shock and along the top of 



the step. Presumably the results of U-PPM and U-PLM 
could be improved by tuning the PPM parameters and 
the 9 parameter in the minmod slope limiter. For this 
problem we use the standard CFL number of 0.5. Reduc- 
ing this number could not improve the results for U-PPM 
and U-PLM. 

A high resolution run employing AMR is shown in 

31 

4.8. Two Dimensional Tests: Shock Tube 

In order to compare with existing multi-dimensional 
SRHD codes we repeat the two-dimensional shock tube 
problem suggest ed bv iDe l Zanna fc Bucci antinil l)2002f) 
and repeated by ILucas-Serrano et alJ l)2004T K This test 
consists of a Cartesian box divided into four equal area 
constant states: 



(p,v x ,v y ,p) NE = 

(p,V x ,Vy,p) NW = 
(p,V X ,Vy,p) SW = 
(p,V X ,Vy,p) SE = 



(0.1,0,0,0.01) 
(0.1,0.99,0,1) 
(0.5,0,0,1), 
(0.1,0,0.99,1) 



where NE stands for the upper right quarter of the box 
(0.5 < x < 1.0, 0.5 < y < 1.0), NW the upper left 
quarter (0.0 < x < 0.5, 0.5 < y < 1.0), SW the lower 
left quarter (0.0 < x < 0.5, 0.0 < y < 0.5) and SE 
the lower right quarter (0.5 < x < 1.0, 0.0 < y < 0.5). 
We use constant zoning of 400 x 400 zones, a T = 5/3 
adiabatic equation of state, outflow boundary conditions 
in all directions and a CFL number of 0.5. 

In the F-PLM run of this test, the parameter 9 in the 
minmod slope limiter is set to 1.2 instead of the default 
value 1.5 to avoid crashes. 

Our results of four schemes are shown in Fig. |H1 
The interface between panels NW and SW and the 
interface between panels SE and SW are stationary 
contact discontinuities with jumps in transverse veloc- 
ity. We note that the contact discontinuities have 
been smeared and numerical artifacts are clearly shown 
in the density contour s at panel SW, as in the re- 
sults in the literature llDel Zanna fc Bucciantinil l2002t 
ILucas-Serrano et al.ll2004[L Indeed less prominent nu- 
merical artifacts also exist in panels NW and SE, though 
they do not show up in the figures of 30 contours. The 
appearance of two curved shocks and the elongated diag- 
onal spike of density in bet ween in panel NE is i n agree- 
ment with the results of ILucas-Serrano et aD (2004), 
whereas the diagonal feature is much less pro minent in 
the results of IDel Zanna fc Bucciantml ((2002) . We note 
that the initial interface between panels NE and NW and 
the interface between panels NE and SE are not simple 
shock waves. A similar test problem with simple waves at 
all initial interfaces (two conta ct discontinuities and two 
shocks ) has been performed bv lMignone. Plewa fc Bodol 
(2005) to test their multidimensional relativistic PPM 
code. 

5. ADAPTIVE MESH REFINEMENT 

Solutions to hyperbolic partial differential equations 
(PDEs) are frequently smooth in large fractions of the 
computational volume yet contain sharp transitions in 
localized regions. In the smooth regions, relatively coarse 
numerical zoning may be sufficient to accurately rep- 
resent the solution, while finer zoning is needed where 
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o.o I U J i i I L J , , I 

0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 

Fig. 8. — Two dimensional shock tube problem at f = 0.4. A 
r = 5/3 adiabatic equation of state was used with outflow bound- 
ary conditions at all boundaries, 400 X 400 Cartesian zones and a 
CFL number of 0.5. Results of F-WENO, F-PLM, U-PPM, and 
U-PLM are shown in panels (a), (b), (c), and (d), respectively. 
Thirty equally spaced contours of the logarithm of proper density 
are plotted. 



sharp transitions occur Adaptive mesh refi nement 
lIBerger fc QgligerlllQSl IBerger fc CoTelTal ll989) creates 
finer numerical meshes to adequately resolve steep gradi- 
ents and thin layers where they exist. For computational 
efficiency, coarser meshes are used in smooth regions so 
that computing power is not wasted over-resolving re- 
gions where high resolution is unnecessary. As the so- 
lution evolves, the mesh structure adapts to provide ap- 
propriate resolution where needed at a given time. 

The SRHD equations are PDEs which admit the for- 
mation of extremely sharp transitions in the form of 
shock waves, thin shells and contact discontinuities. 
These structures coexist with smoothly flowing regions 
in different parts of the computational domain. AMR 
is a powerful tool to handle the resolution required to 
adequately capture these features. 

In the RAM code, we utilize the AMR to ols in the 
FLASH code version 2.3 <|Frvxell et all 12000). which in 
turn is a modified version of the PARAMESH AMR 
package l|MacNeice et al.ll2000j) . PARAMESH is a block- 
structured AMR package. It uses a hierarchy of nested, 
logically Cartesian blocks which typically have eight 
zones per dimension for a total of 8 d zones per block, 
where d = 1, 2 or 3 is the dimensionality of the simu- 
lation. Finer level blocks are a factor of two higher in 
resolution in each direction so that each block is either 
at the highest level of refinement or contains 2 d daughter 
blocks. Flux conservation at jumps of refinement is im- 
posed by replacing fluxes computed at the courser level 
of refinement with appropriate sums of fluxes at the finer 
level. 

Refinement or derefinement of a block is generally de- 
termined by calculating an approximate numerical sec- 
ond derivative of fluid variables which can be specified 
at runtime. Other kinds of criteria can also be used 
for specific problems (see e.g., § 16.41) . In FLASH, the 
one-dimensional normalized second derivative, which is 



a measurement of error, is given by 2 , 

E _ \u i+ 2 - 2Uj + Uj- 2 \ 

\u i+2 - Ui\ + \lti - tti_ 2 | + e(K+2l + 2 KI + K-2I) ' 

(30) 

where the last term in the denominator is a low pass fil- 
ter to avoid excessive refinement on small fluctuations, 
and e is an adjustable parameter for the low-pass filter. 
By defa ult, we use e = 0.01 except for the test problem 
in § 16.11 The one-dimensi onal expression can be general- 
ized to multi-dimension l|Frvxell et alJl2fK)fl) . Typically 
pressure, density and Lorentz factor are used to estimate 
the error norm, E, once every two steps. If the maxi- 
mum error norm on a block is larger than the value of a 
parameter, E re { , the block will be marked for refinement. 
If the maximum error norm on a block is less than the 
value of a parameter, E^eref, the block will be marked 
for derefinement. The default values for the above two 
parameters in the tests shown in this paper are set to 
E le f = 0.8, and Sdcrcf = 0.2 except for the test problem 
in § 16.11 When AMR refines a block, quadratic inter- 
polation algorithms are used for prolongation operations 
between refinement levels. The prolongation operation 
is performed on the conserved variables. Then the con- 
version of conserved variables to primitive variables is 
performed. The conversion could fail in principle, but 
this has never happened in our simulations because the 
refinement is performed when the solution is still very 
smooth. If the conversion fails, the problem can be fixed 
by using first-order prolongation. The solutions on par- 
ent blocks, which are not evolved, are obtained by re- 
striction operations. 

In § 13.31 we discussed the fall-back mechanism, which 
makes our code more robust. When AMR is used, the 
cells in which unphysical results appear are almost al- 
ways at the finest level of refinement. Thus we only need 
to apply the more diffusive schemes for reconstruction 
to the finest level. For coarser levels, we still use the 
standard reconstruction scheme in the recalculation. 

FLASH handles parallelization using the message pass- 
ing interface (MPI) library and uses an estimate of the 
work per processor to balance the computational load 
among processors. 

6. TESTS WITH ADAPTIVE MESH REFINEMENT 

In this section, we present results of the AMR on some 
test problems. All the tests shown are performed with 
the F-WENO scheme, which is very robust and accurate 
for all the tests we have performed. We use F-WENO-A 
and F-WENO-U to denote the F-WENO scheme with 
adaptive mesh and uniform grid, respectively. 

6.1. One-Dimensional Riemann Problem With 
Transverse Velocity: Hard Test 

In § 14.41 we tested out schemes on a one-dimensional 
Riemann problem with non-zero transverse velocity. 
That test problem is relatively easy and can be resolved 
with modest resolution (e.g., 400 zones). To fully ex- 
ercise the code, we perform a very severe test requiring 
very high resolution to resolve the complicate structure 
of the transverse velocity. 

2 Note that the expression for the error norm in FLASH 2.3, 
on which our RAM code is based, is slightly different from that in 
IFrvxell et aT] pQOOT) 
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In this test, the one-dimensional numerical region 
(0 < x < 1) initially consists of two constant states: 
PL = 1000.0, PL = 1.0, v xL = 0.0, v yL = 0.9 and 
Pit = 10~ 2 , pr = 1.0, v xR = 0.0, VyR = 0.9, where L 
stands for the left state, and R the right state. The fluid 
is assumed to be an ideal gas with an adiabatic index 
T = 5/3. The initial discontinuity is at x = 0.5. The 
results at t — 0.6 are shown in Fig. In this test, the 
breakup of the initial discontinuity evolves into a shock 
moving towards the right, a rarefaction wave moving to- 
wards the left, and a contact discontinuity in between. It 
is interesting to note the differences among this test prob- 
lem, a similar problem in § 14.21 which has no transverse 
velocity, and the problem in § 14.41 which has transverse 
velocity on the right side. The presence of large trans- 
verse velocity on the left make the shock move slower 
than that in the previous problems. Compared to the 
previous problems, the density jump is smaller and the 
post shock dense shell is wider. However, these differ- 
ences do not make this problem easier. The distance 
between the tail of the rarefaction wave and the contact 
discontinuity is much smaller in this test. More inter- 
esting is the structure of the transverse velocity. Inside 
the rarefaction wave, the transverse velocity increases 
from v y = 0.9 at the head of the rarefaction wave to 
v y = 0.9602 and then decreases to v y = 0.9472 at the 
tail of the rarefaction wave. While the Lorentz factor 
increases monotonically from W = 2.29 at the head to 
W = 35.8 at the tail of the rarefaction wave. Both the 
normal and transverse velocity stay constant in the thin 
shell between the tail of the rarefaction wave and the con- 
tact discontinuity. The presence of a very large Lorentz 
factor (~ 36) in that thin shell on the left side of the con- 
tact discontinuity requires extremely high resolution to 
resolve the structure. Across the contact discontinuity, 
pressure and normal velocity do not change, while den- 
sity and transverse velocity have a jump in their values. 
At the shock front, the post shock state has a transverse 
velocity of v y = 0.7721, which is smaller than the pre- 
shock value v y = 0.9. 

The results (Fig. EJ) clearly demonstrate the demand of 
high resolution when strong shear flows present. In low 
resolution runs, the position of the shock front and the 
contact discontinuity were not correctly captured. Even 
with 51200 zones, there is still visible errors for trans- 
verse velocity at the contact discontinuity. In the cases 
like this, we argue that AMR could be a powerful tool. 
For the sake of comparison, we run the test with vari- 
ous uniform grids and AMR grids with correspondingly 
equivalent resolutions. For this test, the AMR parame- 
ters (0 are set to e = 0.005, E le f = 0.5, and -Edcrcf = 0.1. 
For the AMR grids, the fine zoning follows the shock 
front and the contact discontinuity automatically, while 
coarse zoning is used at smooth region like the middle 
of the rarefaction wave. In one of the runs, 8 levels of 
refinement are used. That is the finest zones are 128 
times smaller than the coarsest zones. The AMR makes 
the resolution of this run equivalent to 51200 uniform 
zones. Meanwhile AMR runs are much faster than their 
corresponding uniform zoning runs. The measurement 
of global errors is shown in Table As we expected, 
results from AMR are comparable to those from uniform 
zoning, even though coarse zoning is used in most of the 
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Fig. 9. — ID Ricmann problem with transverse velocity on both 
sides at t = 0.6. Results of F-WENO-A with three different res- 
olutions arc shown. The lowest refinement level is 1 in all three 
runs, while the highest levels are 1, 4, and 8, for top, middle, and 
bottom panels, respectively. The level 1 contains 400 zones. Thus 
the equivalent numbers of zones are 400, 3200, and 51200, for top, 
middle, and bottom panels, respectively. Color lines show the nu- 
merical results for proper density (red), pressure (green), normal 
velocity (blue), and transverse velocity (cyan), whereas solid lines 
show the exact solutions. 



TABLE 7 

Li ERRORS OF THE DENSITY FOR THE ID RlEMANN PROBLEM 
WITH TRANSVERSE VELOCITY ON BOTH SIDES. 



Mesh 


Levels a 


N h 


L\ Error 


Convergence Rate 


Uniform 


1 


400 


5.21e-l 






1 


800 


3.63e-l 


0.52 




1 


1600 


2.33e-l 


0.64 




1 


3200 


1.26e-l 


0.89 




1 


6400 


6.49e-2 


0.96 




1 


12800 


3.38e-2 


0.94 




1 


25600 


1.80e-2 


0.91 




1 


51200 


9.95e-3 


0.86 


Adaptive 


1 


400 


5.21c-l 






2 


800 


3.63e-l 


0.52 




3 


1600 


2.33e-l 


0.64 




i 


3200 


1.26e-l 


0.89 




5 


6400 


6.55e-2 


0.94 




6 


12800 


3.49e-2 


0.91 




7 


25600 


1.90e-2 


0.88 




8 


51200 


1.07c-2 


0.83 



a Total levels of the grid 

b For adaptive meshes, this means the equivalent number of 
zones. 



computational region. This is due to the fact that the 
global error is dominated by the discontinuities where 
high resolution is used in AMR runs. 

Note that the behavior of transverse velocity in one- 
dimensional Riemann problems is purely due to relativis- 
tic effects. One might have the following argument which 
would lead to the wrong conclusion that transverse veloc- 
ity cannot change in one-dimensional Riemann problems. 
The initial transverse velocity before the breakup could 
disappear by transferring to a new reference frame which 
moves along the transverse direction in the original lab- 
oratory frame. In the new frame, no transverse velocity 
would develop if the initial transverse velocity in that 
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frame is zero. After the breakup of the initial disconti- 
nuity, one could transfer back to the original frame. The 
second Lorentz transformation would make the trans- 
verse velocity in the original back to its original value. 
Thus the transverse velocity in the original frame should 
be fixed during the evolution of the Riemann problem. 
However, relativistic effects invalidate the above argu- 
ment. Suppose the original setup is a two-dimensional 
tube and there is a diaphragm at the initial discontinu- 
ity. The breakup of the initial discontinuity is caused by 
taking the diaphragm away. This problem can be con- 
sidered one-dimensional in the original frame because of 
the obvious symmetry. However, in the second reference 
frame, which moves along the transverse direction, the 
problem can no longer be considered one-dimensional. 
The simultaneous disappearance of the diaphragm at 
different places in the transverse direction in the orig- 
inal frame does not happen simultaneously in the sec- 
ond frame. Therefore, transverse velocity does change in 
both the original and the second frame. This conclusion 
can also be seen directly from the governing equations 
(Eq. |3J and the definition of the conserved mass and mo- 
mentum densities (Eqs. 0and|H}. These imply that the 
SRHD evolution equations preserve hWvt across shocks 
and rarefaction fans, where Vt is a velocity component 
transverse to the propagation direction of these waves. 
Since, in general, h and W change, Vt is generally not 
conserved. 

6.2. Advection Across Jumps in Refinement 

On AMR grids, there are jumps in resolution at the 
boundaries between computational regions at different 
levels of refinement. Spurious effects such as the reflec- 
tion of waves could arise at these boundaries. If AMR 
performs as desired, the effects should be negligible be- 
cause the structure at the refinement boundaries should 
be very smooth. To assess the effects of jumps in refine- 
ment in the worst case scenario, we present an advection 
test on a fixed staggered mesh. We set up a compu- 
tational grid with a range of refinement levels, but we 
turn off adaptive refinement and derefinement so that 
the mesh remains spatially variable but constant in time 
(see Fig. EU). 

In this test, the computational region consists of a two- 
dimensional box in Cartesian coordinates with —0.45 < 
x < 0.45 and —0.45 < y < 0.45. All the boundaries are 
periodic. A fixed mesh with 4 levels of refinement is set 
up on the grid. On this non-adaptive staggered mesh 
(Fig. ITU)), there is a reference state of p — 1.0, p lc f = 1.0, 
v x = 0.72, v y = 0.54, with a pulse only in density. The 
density profile is given by p (r) (Eqs.ED&EJ, where 
r is the distance to the center of the box, and the two 
parameters (the amplitude and width of the pulse) in 
the density profile are set to a — 10.0, L — 0.2. The 
adiabatic index in the equation of state is T = 5/3 in 
this test. The density pulse will advect on the periodic 
numerical grid at an angle of arctan(4/3) with respect 
to the x-axis and move across the refinement boundaries. 
The results are shown in Fig. As it propagates, the 
pulse becomes wider due to numerical viscosity. How- 
ever, no spurious waves are visible, and the pulse is still 
symmetric. At t = 10, the fluctuation in pressure is only 
at the level of 10 -13 . Such a small fluctuation, which is 
probably due to round-off errors in floating point opera- 
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Fig. 10. — Advection across jumps in refinement. F-WENO 
with 4 levels of refinement fixed in time is used in this calculation. 
The block structure of the mesh is shown by the black solid lines. 
Each block contains 8x8 zones. Results of the density structure at 
t = 0.0 (upper left panel), t = 3.5 (upper right panel) t = 7.5 (lower 
left panel), and t = 10.0 (lower right panel) are shown. The density 
pulse advects at a speed of v = 0.9 at an angle of arctan (4/3) with 
respect to the ir-axis and moves across the refinement boundaries.. 
The boundaries of the numerical box are periodic. 



tions, suggests that spurious effects due to the jump in 
refinement is acceptably small in our simulations. 

6.3. Relativistic Jet In Two-Dimensional Cylindrical 
Geometry 

Many authors have performed two-dimensional sim- 
ulations to study the morphology and dynamics of 
relativistic 
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AGNs. These results confirmed the formation of the 
basic features of a jet, i.e., beam, cocoon, Mach disk 
and bow shock, which have bee n observed in Newton ian 
hydrodynamic simulations fe.g.. lNorman et aDll982|) . It 
was also found that relativistic jets are more stable and 
propagate more efficiently than Newtonian ones. Three- 
dimensional relativistic hydrodynamic simulations of 
jets have also been performed to study three-dimensional 
effects, such as 3D instabilities and precession of jets 
(Hughes eTaLll200l lAlov et alJl2003^ . 

In this test, we simulate a relativistic jet which is rele- 
vant in astrophysics. The computational region is a two- 
dimensional cylindrical box (0 < r < 15, < z < 45). 
The details of treating curvilinear coordinates including 
both cylindrical and spherical coordinates can be found 
in appendix For the sake of comparison, we use 
the same paramet ers for this problem as model C2 of 
iMarti et al.l 1^997). The initial parameters of the rela- 
tivistic jet beam are, Vb = 0.99 and pb = 0.01. The classi- 
cal Mach number of the jet is set to Mb — 6, and the cor- 
responding relativistic Mach number, Mb = WbMb/W s , 
where Wb is the Lorentz factor of the jet beam and 
W s is the Lorentz factor of sound speed in the jet, is 
~ 42. r = 5/3 is used for the equation of state. Ini- 
tially the computational region is filled with a uniform 
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medium with p m = 1, and v m — 0. The pressure of the 
medium, which is chosen to match the pressure of the 
jet, is p m — pb = 0.000170305. The initial relativistic jet 
is injected through the inlet part (r < 1) of the low-z 
boundary by assigning the state of jet material to the 
boundary. Outflow boundary conditions with zero gra- 
dients of variables are used at the other part (r > 1) of 
the low-z boundary, high-z boundary, and high-r bound- 
ary. A reflecting boundary is used at the symmetric axis 
where r = 0. 

Our numerical results at t = 100 are shown in Fig. 1111 
It is shown that no carbuncl e, which is a pathological 
phenomenon in some schemes (|Quirklll994|) , is generated 
in these simulations. We have performed the simula- 
tions with two different resolutions using F-WENO-A. 
In both calculations, the lowest level of the grid consists 
of 24 x 72 zones. The corresponding spatial resolution 
at this level is ~ 1.6 zones per jet beam radius. The 
total levels of refinement are 5 and 7, for the two calcu- 
lations, respectively This corresponds to an equivalent 
resolution of ~ 26 and ~ 102 zones per jet beam ra- 
dius, for the two calculations, respectively. The model 
in this test is highly supersonic, and relativistic effects 
from ultra-rel ativistic motion d ominates those from in- 
ternal energy ijMartf et all! 9971. The expected morpho- 
logical features of such a relativistic jet are observed. A 
bow shock is formed due to the supersonic motion of the 
jet. The medium is shocked by the bow shock. The jet 
beam is slowed down at the Mach disk and feeds the co- 
coon. The shocked jet beam moves sideways and then 
even backwards. The discontinuity between the shocked 
jet material and shocked medium material admit Kelvin- 
Helmholtz instabilities in the cocoon. The average speed 
of the jet head is ~ 39. This agrees wit h the essen- 
tially one-dimensional analytic estimate, 42 ijMarti et alJ 
1997). Comparing the results from the two calculations 
with different resolutions, we found that the global struc- 
ture and the propagation speed of the jet are almost 
identical though more mixing due to Kelvin-Helmholtz 
instabilities is observed in the high resolution run. 

The presence of the strong shear motion between the 
jet beam and the back flows challenges numerical codes. 
Even though a smaller CFL number, 0.3, is used in this 
test, our code still falls back to more diffusive schemes 
occasionally. We consider this as a signal that the cal- 
culations are still under-resolved for such strong shear 
flows. We note that this is the only test that the fall- 
back scheme (§ 13. 3|) is used. 

Admittedly, this test problem is not ideal for AMR. At 
later times, most of the region is at the highest refinement 
level because of the rich structures inside the bow shock. 
However, AMR still saves a lot of computer time because 
the external medium requires very low resolution during 
most of the time. 

6.4. Relativistic Kelvin-Helmholtz Instability 

The relativistic Kelvin-Helmholtz instability is of great 
interest in many astrophysical problems. For instance, 
a key unanswered question for GRB jets is whether or 
not the relativistic outflows can remain sufficiently clean 
or whether baryons will be mixed into the flow lower- 
ing the maximum asymptotic Lorentz factor below the 
values required by observations. This problem is espe- 
cially relevant for collapsars in which the jet must prop- 
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Fig. 11. — Highly supersonic relativistic jet in 2D cylindrical 
geometry at t = 100. F-WENO-A is used in the two calculations 
with (a) 5 and (b) 7 levels of refinement. The equivalent resolutions 
are (a) 384 X 1152 and (b) 1536 X 4608 zones. This corresponds to 
a resolution of ~ 26 and ~ 102 zones per jet beam radius, for (a) 
and (b), respectively. The CFL number used in this test is 0.3. 



agate through a star. The Kelvin-Helmholtz instabil- 
ity for relativistic jets has been studied both numeri- 
cally and analyt i cally (e.g., iRosen et~a II IT999T: IHarded 
1200(1 I.Xlov et all 12002). More recently, iPerucho et al 
H2004albL 12005^ have performed a series of studies on 
the linear growth and long-term nonlinear evolution of 
two-dimensional relativistic planar jets. In this section 
we present two numerical tests of the Kelvin-Helmholtz 
instability. 

It is well known that ultra-relativis tic flows suppress 
the Kelvin-Helmholtz instabilities (see lBodo et alJl2004, 
for a recent analytical work). In the first test in this 
section, we present a numerical simulation of the Kelvin- 
Helmholtz instability with mildly relativistic speeds to 
demonstrate the ability of AMR for problems involving 
small scale structures. The computational region con- 
sists of a two-dimensional box with < x < 1 and 
^5 < y < 5. F-WENO-A with 5 levels of refinement is 
used in this calculation. The equivalent spatial resolution 
is 1024 x 10240 zones. The upper half of the box is filled 
with a gas with p = 1, p = 1000, v x = 0.9 and v y = 0, 
whereas the bottom half with p = 10, p = 1000, v x = 
and v y = 0. V = 5/3 is used for the equation of state of 
ideal gas. Periodic boundary conditions are used for the 
right and left boundary. The initial interface which sep- 
arates the two fluids is described by y — 0.01sin(27ra;). 
This small perturbation triggers the Kelvin-Helmholtz 
instability. Fig. ^| shows a series of snapshots of the 
results. The rolling up of the interface forms small vor- 
tices. This process continues to form larger vortices from 
smaller vortices and becomes unstable. 

To allow for a direct comparison with previous re- 
sults on the relativistic Kelvin-Helmho l tz inst ability, we 
repeated model D10 in IPerucho et alJ {2005) for a rel- 
ativistic sheared planar jet using the RAM code. In 
this test, the computational region consists of a two- 
dimensional box in Cartesian coordinates with < x < 8 
and — 40 < y < 40. The jet and ambient medium initially 
have the same pressure, po — 2.0. The jet, which moves 
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Fig. 12. — Rclativistic Kclvin-Helmholtz Instability. Results at 
t = 0, 1, 2, 3, 4, and 5 are shown in panels (a), (b), (c), (d). 
(e), and (f), respectively. F-WENO-A with 5 levels of refinement 
is used in this calculation. The equivalent spatial resolution is 
1024 X 10240 zones. The whole computational region consists of 
a two-dimensional box with < x < 1 and —5 < y < 5. In this 
figure, only part of the domain is shown. Initially, the upper half 
of the fluid moves at a speed of 0.9c, while the bottom half is at 
rest. 



along the x-direction, has initial density of pj =0.1, and 
an initial Lorentz factor of Wj — 10.0. The static ambi- 
ent medium has initial density of p a = 1.0. To allow for 
a continuous transition between the jet and medium, the 
profiles of density and velocity are smoothed as follows, 



po(y) = Pa 



Pa - P] 



cosh (y/RjY 



vx{y) 



cosh (y/Rj 



(31) 



(32) 




where the jet radius Rj is set to 1.0, and the steepness 
parameter m is set to 25. The adiabatic index in the 
equation of state is T = 4/3. A periodic boundary con- 
dition is imposed at both the low- a; and high- a; boundary. 
To excite the Kelvin- Helmholtz instability, a small initial 
perturbation is given to the transverse velocity, v y : 

N-i ~\ 
^ sin [(n + l)k n x + </>„] sin 2 [(n + l)wy] y- ? > 

V/-1 ~j 

sin [(m + l)k m x + <f> m ] sin 2 [(m + l)iry] > , 

m=0 ) 

(33) 

where v\ = 5.77 x 10 -6 (Perucho 2005, private commu- 
nication) is the amplitude, fc m , n are the wavenumbers of 
the modes, 4> m ,n are the random phases for each mode, 
N = 4 is the number of symmetric modes, and M = 4 is 
the number of antisymmetric modes. The wavenumbers 
for these modes are given by kj = (j + l)2n/L, where 
j = 0, 1,2,3, and L = 8.0 is the length of the periodic 
box in x-direction. 

In the test of model D10, 6 levels of mesh refinement 
are used. On the finest level, the resolution is 32 zones 
and 256 zones per unit length at x and y-direction re- 
spectively. In this test, resolution is concentrated at the 
interface between the jet material and ambient medium. 
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Fig. 13. — Evolution of the relative amplitude of perturbations 
for model D10 of relativistic Kelvin-Hclmholtz instability. The am- 
plitudes of perturbation in the "jet reference frame" are defines as 
0.5(1)3;, max -fa;,min)i and 0.5(i>j,, m ax — v v , mm ) , for the longitudinal 
velocity (solid line) and perpendicular velocity (dot line) respec- 
tively. The amplitude of the pressure perturbation (dash line) is 
defined as (p max — Po)/pOi where po is the initial pressure. 



Instead of using the general refinement criteria (§0), we 
use the composition of jet material and ambient medium 
to determine whether refinement or derefinement should 
be performed. More specifically, the refinement level will 
be brought to the highest if the mass fraction of jet ma- 
terial is 0.001 < Xj < 0.999, or otherwise will be marked 
for derefinement. 

FigurcslT3lfell4lshow the results of mod el DIP using F- 
WEN O-A. Our results agree with those of lPerucho et alJ 
(2005). The evolution of the Kelvin-Hclmholtz instabil- 
ity consist s of three regimes: li near, saturation and non- 
linear (see iPerucho et aT1l2005l for more details). In the 
linear regime, the perturbations of the longitudinal ve- 
locity, v x , transverse velocity, v y , and pressure all grow 
exponentially (Fig. 113(1 . At the end of the linear regime, 
tn n w 235, the growth of the longitudinal velocity per- 
turbation starts to saturate, while both the transverse 
velocity and pressure perturbations still grow exponen- 
tially. In the saturation regime, the continuous growth of 
the perturbation of transverse velocity, v y results in the 
distortion of the jet beam (Fig.ll4|) . The transition of the 
saturation to fully non-linear happens when the pressure 
perturbation and transverse velocity perturbation reach 
their peak at i p0 ak ~ 340 and t sat ~ 355 respectively. 
In the non-linear regime the jet material is completely 
mixed with the ambient medium by turbulent motions. 

6.5. Emery Step 

The Emery step is a standard test involving the flow in 
a wind tunnel encountering a reflecting step boundary. 
We have shown the results with four different schemes (F- 
WENO, F-PLM, U-PPM, and U-PLM) on uniform grids 
in § 14.71 Here, we repeat the test with F-WENO-A on 
adaptive meshes. The setup of this problem is the same 
as in § 14.71 except that an AMR grid with 5 levels is used 
instead of uniform grids. The equivalent resolution is 
3840 x 1280 zones on a two-dimensional box of < x < 3 
and < y < 1. In other words, this corresponds to a 
zoning of Ax — Ay — 1/1280. Again, no special treat- 
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Fig. 14. — Density structure of model D10 from Pcrucho ct al. 
120051) of the relativistic Kelvin-Helmholtz instability. F-WENO-A 
with 6 levels of refinement is used in this calculation. The whole 
computational region consists of a two-dimensional box, < x < 
8, —40 < y < 40, with periodic boundary conditions in the re- 
direction. Results at t = 100, 235, 355, and 600 are shown. The 
four panels show the proper density structure during the middle 
of the linear regime (upper left), at the end of the linear regime 
(upper right), at the end of the saturation regime (lower left), and 
during the fully non-linear regime (lower right), respectively. The 
solid white lines denote where the mass fraction of the jet material 
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merit for the step other than a reflecting boundary is 
used because we are only interested in the global struc- 
tures of shocks. This problem involves a Mach reflec- 
tion of a strong shock by the upper boundary, and then 
the reflected shock is reflected again by the step. AMR 
should work very well for this problem because much of 
the computational volume is smooth flow requiring only 
low resolution. Thus this problem can test not only the 
ability of the AMR code to capture strong shocks, but 
also its ability to adaptively refine and derefine the mesh 
as discontinuities evolve. It is shown in Fig. ^| that the 
AMR in our code captures the sharp transitions in the 
flow where and when they occur. Even the weak contact 
discontinuity originating from the bottom of the Mach 
stem is detected and captured by the AMR. It is inter- 
esting to note that no Kelvin-Helmholtz instability de- 
velops along the contact discontinuity though it does for 
non-relativistic flow. 

6.6. Double Mach Reflection Of A Relativistic Shock 

Double Mach reflection of a strong shock has often 
been used to test Newtoni an hydrodynamics cod e s sinc e 
it was first introduced bv iWoodward fc Colellal |l984). 
The computational region consists of a two-dimensional 
box with < x < 4 and < y < 1. A reflecting wall 
is placed at the x > 1/6 part of the low-y boundary. A 
strong shock is moving towards right with its front mak- 
ing a 60 degree angle with respect to the x-axis. The 
reflection of the incident shock could give arise to a dou- 
ble Mach reflection. Two Mach stems and two contact 
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Fig. 15. — Emery step at t = 4.0 using F-WENO-A with 5 levels 
of refinement. The equivalent resolution is Ax = Ay = 1/1280. 
Proper mass density and pressure are plotted in the upper and 
bottom panel, respectively. The velocity field is also plotted in 
the upper panel. In the bottom panel, the block structure of the 
adaptive mesh is shown. Each block contains 8x8 zones. 



discontinuities would appear. The first contact discon- 
tinuity extends from the first triple shock point to the 
reflecting wall and then is pushed by the pressure gradi- 
ent there and forms a jet parallel to the wall. 

The original setup of the problem cannot be used for 
special relativistic hydrodynamics because the specific 
parameters chosen are not appropriate for the relativis- 
tic case. To test our relativistic schemes, we modified 
the original parameters of the problem. The adiabatic 
index in the equation of state is T = 1.4. The unshocked 
gas has po — 1.4, po = 0.0025, and vq = 0. The classical 
Mach number of the shock, M — v s /c s , where v s is the 
speed of the shock front and c s is the sound speed of the 
unshocked gas, is 10. Using the relativistic shock jump 
conditions, we can set up the shocked gas. The shock 
front initially intersects the lower boundary at the edge 
of the wall, x = 1/6. The x < 1/6 part of the low-y 
boundary and the low- a; boundary are set to the exact 
post-shock state. At y = 1, the boundary conditions are 
set to either post-shock or pre-shock conditions depend- 
ing upon the exact motion of the shock. 

The post-shock gas moves at a speed of v\ m 0.4247, 
and the shock speed is v s ~ 0.4984. We would have 
made the shock ultra-relativistic. However, an ultra- 
relativistic shock cannot generate a Mach reflection. This 
is qualitatively understandable. Suppose that the shock 
front makes an angle of 9 with respect to the wall and 
the shock velocity is v s . The velocity of the intersection 
point of the shock and the wall would be v s /sm(9). With 
ultra-relativistic shock velocity, the intersection point 
can move faster than the speed of light. Note that this 
does not violate relativity because the velocity of the 
intersection is not a physical velocity. Even if a Mach 
stem can be formed initially, the vertical Mach stem can- 
not move faster then the speed of light, and thus cannot 
follow the motion of the oblique incident shock. There- 
fore no permanent Mach stem is possible when the shock 
moves too fast. 
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Fig. 16. — Double Mach reflection of a relativistic shock at t = 4. 
Results with and without AMR are shown in panels (a) and (b), 
respectively. Thirty equally spaced contours of proper density are 
plotted. Only part of the computational domain is shown. Results 
from (a) F-WENO-A and (b) F-WENO-U are almost identical. 



We have performed this test with both adaptive mesh 
(F-WENO-A) and uniform grid (F-WENO-U). Our re- 
sults are shown in Fig.^j] For the adaptive mesh, 4 levels 
of refinement are used with the lowest level containing 
64 x 16 zones. The equivalent resolution is 512 x 128. 
For the uniform grid, 512 x 128 zones are used. The 
results from F-WENO-A can be compared with those 
from F-WENO-U. The AMR scheme works very well 
and the highest level of refinement follows the motion 
of shocks and discontinuities while most of the compu- 
tational region is at lower refinement levels. It is shown 
in Fig. El that density contours for F-WENO-A and F- 
WENO-U are almost identical. Our results do not suf- 
fer from some pathological behaviors obtained wi th some 
schemes (|Woodward fc Colellall98llQuirkll99^ . There 
are no kinked Mach stems. And the region behind the 
slowly moving shock near the left edge of the reflecting 
wall is very smooth. 

6.7. Three-Dimensional Blast Wave With Spherical 
Symmetry 

In this test, we study a three dimensional spherical 
blast wave using Cartesian coordinates. Since no ana- 
lytic solution exists for this problem, we compare the 
three dimensional solution to a high resolution one di- 
mensional simulation run with spherical coordinates. For 
the sake of comparison, we perfo rm the same test run by 
iDel Zanna fc Bucciantml <|2002ft . See their Fig. 8. 

The computational domain is a cubic box, < x < 1, 
0<y<l, < z < 1, with a base resolution of 40 zones 
in each direction. We use up to 4 levels of AMR so the 
effective resolution is 320 x 320 x 320 zones. An initial 
discontinuity is at R = 0.4, where R is the distance to 
the center of the Cartesian coordinates. The region in- 
side the interface, R < 0.4, contains a gas with p = 1, 
p = 1000, whereas the region outside it, R > 0.4 con- 
tains a gas with p = 1, p = 1. r = 5/3is used for the 
equation of state. Both gases are at rest initially. Simi- 
lar to the one-dimensional Riemann problem, the decay 
of the initial discontinuity will give arise to a spherically 
symmetric shock moving outwards, a spherical rarefac- 
tion wave moving inwards, and a "contact discontinu- 
ity" in between. Note that the system could eventually 
evolve into a selfsimilar Blandford-McKee solution for 
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Fig. 17.— Three-dimensional blast wave at t = 0.4. F-WENO-A 
with 4 levels of refinement is used for the simulation. The equiva- 
lent resolution is 320 X 320 X 320 zones for a cubic box, < x < 1, 
< y < 1, < 2 < 1. The results along the major diagonal of 
the cubic box are plotted. Three-dimensional numerical results are 
shown in symbols, whereas the high resolution numerical results in 
one-dimensional spherical coordinates are shown in solid lines. We 
show proper mass density (square), pressure (triangle) and total 
velocity (plus sign). 



ultra-relativistic blast waves JBlandford fc McKedllTm!) 
if the initial conditions are modified to make an even 
stronger explosion. The boundaries at x = 0, y = and 
z = are reflecting, whereas all other boundaries are 
zero gradient outflows. 

Fig. El shows the result at t = 0.4 for a simula- 
tion using a CFL number of 0.3. The solid line is a 
one-dimensional simulation run using spherical coordi- 
nates with 4000 uniform zones. The details of the treat- 
ment of curvilinear coordinates can be found in the ap- 
pendix. The results with three-dimensional Cartesian 
coordinates agree with the high resolution results with 
one-dimensional spherical coordinates. The shock front 
is captured with ~ 3 zones. The spherical symmetry 
is kept in the overall structure of the numerical results 
using Cartesian coordinates. 

6.8. Scaling Tests on Parallel Machines 

Many challenging problems in numerical astrophysics 
require the computing power made possible by massively 
parallel supercomputers. These include state-of-the-art 
machines available through the national supercomputing 
centers and the national laboratories, as well as the Linux 
clusters available in an increasing number of institutions. 
In order to take full advantage of available computational 
power, it is necessary to develop computer code and uti- 
lize packages which run efficiently on a variety of parallel 
platforms and scale well with the number of processors 
used. Large simulations, especially in 3D, require the use 
of large numbers of parallel processors to run in a reason- 
able amount of time. In order to use parallel resources 
efficiently, it is necessary to test how well a code scales 
with number of processors used. Scaling information is 
important for planning the size of numerical simulations 
which are possible to run efficiently. It is also impoor- 
tant for using an efficient number of processors for a given 
job. We have tested RAM on several major national su- 
percomputers as well as smaller clusters and find that it 
scales very well (Figure IT%|) . 

RAM has been extensively tested on NASA's Columbia 
supercomputer on up to 504 processors. A scaling anal- 
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Fig. 18. — Scaling of two "constant total work" parallel sim- 
ulations run on NASA's Columbia supercomputer from 32 to 340 
processors. The top pair of lines are for a simulation with a spa- 
tially variable but temporally fixed mesh. The bottom pair is for a 
simulation using adaptive refinement. The dashed lines represent 
perfect scaling. 



ysis for a stellar collapse simulation run with a spatially 
variable temporally fixed grid shows how RAM scales 
from 32 to 340 processors (Figure IT51 upper lines). For 
this constant total work job, RAM runs at 93% max- 
imum efficiency when going from 32 to 128 processors 
and 88% maximum efficiency when going from 128 to 
340 processors. We define efficiency when going from 
m to n processors as e = t m /t n * (m/n), where t x is 
the run time with x processors. In a different simula- 
tion of a relativistic shock using a spatially and tempo- 
rally adaptive mesh (Figure 1181 lower lines) , processing 
time scales with 80% of maximum efficiency when going 
from 32 to 128 processors, and 65% efficiency going from 
128 to 340 processors on Columbia. This simulation re- 
quired more communication among processors due to the 
continual refinement and derefinement as flow features 
moved through the computational domain. 

Since the RAM code is based on FLASH 2.3, we expect 
the scaling and efficiency to be similar. More details 
on the s caling and efficiency of the FLASH code can be 
found in lFrvxell etaD l|2000jK 

7. SUMMARY AND DISCUSSION 

We have presented a new special relativistic hydro- 
dynamics code with adaptive mesh refinement. The 
code is modular and includes four combinations of re- 
construction and hydrodynamics solvers we have termed 
F-WENO, F-PLM, U-PPM and U-PLM. In this paper, 
we focus attention on F-WENO, a characteristic- wise, 
finite difference WENO scheme utilizing the full char- 
acteristic decomposition of the SRHD equations. The 
scheme is fifth-order accurate in smooth regions. This 
is the first time that this high-order scheme has been 
implemented for relativistic hydrodynamics. We demon- 
strate that, while somewhat more complex, this method 
is highly accurate and stable. It has the added advantage 
of not containing tunable parameters which substantially 
modify the performance of the algorithm. 

The U-PPM scheme utilizes an approximate Riemann 
solver (e.g., modified Marquina flux formula) an d is 
equivalent to the GENESIS code ijAlov et al]ll999fl . U- 
PLM and F-PLM, use linear reconstruction and are more 
diffusive, though somewhat stabler, versions of U-PPM 



and F-WENO. For some problems these more diffusive 
schemes perform better when suppression of spurious os- 
cillations near discontinuities is desirable. The higher- 
order methods sometimes produce numerical oscillations 
near strong discontinuities leading to values of conserved 
variables which are inconsistent with primitive fluid vari- 
ables. This inconsistency can lead to code crashes. 
We therefore have implemented a failsafe scheme, which 
saves the previous solution at all times. If a numerical 
step results in an unacceptable solution, we return to the 
previous step and repeat the hydrodynamics calculation 
with increasingly diffusive schemes which, due to their 
diffusive nature, are more likely to produce consistent 
variables. This procedure greatly increases code robust- 
ness. 

Numerical relativistic hydrodynamics is more difficult 
than the Newtonian case for several fundamental reasons. 
First, due to relativistic effects, structures in relativistic 
flows tend to be thinner requiring higher spatial resolu- 
tion. Second, unlike shocks in Newtonian hydrodynam- 
ics, the density jump across relativistic shocks can be 
arbitrarily large and is limited only by the Lorentz fac- 
tor (see Eq.EEJl- Third, the non-linear coupling between 
the velocity components in relativistic flows remains a 
very challenging numerical problem. The Lorentz factor 
depends on all velocity components and couples trans- 
verse velocity into the dynamics in the normal direction. 
The difficulty is in resolving contact discontinuities which 
move with respect to the numerical grid. Some smearing 
of contact discontinuities is inevitable due to numerical 
diffusion error. When transverse velocity components 
are present they experience jumps across contact discon- 
tinuities with corresponding jumps in the Lorentz factor, 
unlike the case with no transverse velocity in which the 
Lorentz factor is continuous across contact discontinu- 
ities. Because the conserved variables depend on differ- 
ent powers of the Lorentz factor, the numerically smeared 
states generate spurious waves from the discontinuity 
which corrupt the solution. An additional, though less 
serious, difficulty is that pressure jumps across the shock 
and inside the rarefaction fan have a steeper dependence 
on normal velocity when transverse velocities are present. 

These problems can be solved by increasing spatial res- 
olution. A major advantage of RAM is its high accuracy 
coupled with AMR for extremely high effective resolu- 
tion. We have found that RAM is able to achieve suffi- 
cient resolution to correctly solve challenging numerical 
problems. Since AMR concentrates the spatial resolu- 
tion where it is needed, sufficient accuracy is achieved 
efficiently, enabling challenging multi-dimensional simu- 
lations to be undertaken. 

The high accuracy achieved with RAM is demon- 
strated to be of critical importance for solving relativis- 
tic flows with strong shear, even in the more challeng- 
ing case when transverse velocity is present in high pres- 
sure states into which rarefaction fans are propagating. 
Under-resolved simulations with many combinations of 
reconstruction and hydrodynamics solvers produce in- 
correct post-shock values and positions as seen in Fig. EI 
We believe this problem is generic to most or all rel- 
ativistic hydrodynamics codes currently in use (see also 
Mienonc. Plewa & Bod ol2005[) . In this work, we succeed 
in obtaining agreement with the analytic solution for a 
Riemann problem with transverse velocity of v ~ 0.9 by 
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using 8 levels of AMR. We are confident that we can use 
RAM to adequately resolve flows with strong shear which 
have not previously been adequately resolved. 

It is worth noting that U-PPM and U-PLM schemes 
do not require the complete characteristic information 
of the Jacobian matrix of the SRHD system. How- 
ever, this information is needed for the characteristic- 
wise schemes (F-WENO and F-PLM) making them 
more difficult to implement. Our numerical experi- 
ments have shown that the extra effort is justified. 
We have also implemented other F-X schemes for com- 
parison including t he fifth-order monoticit y-preserving 
scheme (MP5) of iSuresh fc Huvnhl I 1997T) and third- 
order WENO (WENQ3) Piang fc Shul ll996L In our nu- 
merical experiments we have found that MP5 generates 
excessive spurious numerical oscillations, while WEN03 
is exceedingly diffusive. The fifth-order WENO scheme 
implemented in this work (see § 13.1(1 is free from these 
liabilities. 

In order to simulate systems of astrophysical interest 
we have included physics modules relevant for the study 
of GRBs, SN and black hole accretion. RAM includes 
physical equations of state, neutrino cooling and nuclear 
physics. It is intended initially to address open theo- 
retical questions in the study of gamma-ray bursts and 
supernova explosions. 
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APPENDIX 

SPECIAL RELATIVISTIC HYDRODYNAMICS IN CURVILINEAR COORDINATES 

When a method which is valid in Cartesian coordinates is extended to curvilinear coordinates care must be taken 
to avoid the introduction of errors related to the mesh geometry. We treat the extension to non-Cartesian coordinates 
carefully to avoid errors especially near coordinate singularities. Here we describe the extension of our code to 
cylindrical and spherical coordinate geometries. 

In cylindrical coordinates (r, 9, z), the governing SRHD equation can be written as, 
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where the subscripts, r, 9 and z stand for radial, azimuthal and axial directions in cylindrical coordinates. All fluid 
variables have the same meaning as in § [21 
The above equations for cylindrical coordinates can be easily discretized into, 
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where the subscripts i, j and k denote the r-, 9- and z-discretization, respectively. The subscripts i ± 1/2, j ± 1/2, 
and k ± 1/2, refer to cell interfaces. ^Ji,j,k is the mean value of the conserved variable at the cell (i,j,k). F are the 
numerical fluxes at the cell int erfac es. The source term, Sij t k, can be calculated according to the right-hand sides of 
Equations [ATI lA"2l lAll lA"4l and lA"5l 
Similarly, the governing equation in spherical coordinates (r, 9, <fi) can be written as, 
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where the subscripts, r, and </> stand for radial, polar and azimuthal directions in spherical coordinates. All fluid 
variables have the same meaning as in § [5] 
The discretized equations for spherical coordinates read 
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where the subscripts i, j and k denote the r-, 0- and (^-discretization, respectively. The subscripts i ± 1/2, j ± 1/2, and 
k ± 1/2, refer to cell interfaces. 
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